Distributed Temperature Measurement in a Self-Burning Coal Waste Pile through a GIS Open Source Desktop Application

Geographical Information Systems (GIS) are often used to assess and monitor the environmental impacts caused by mining activities. The aim of this work was to develop a new application to produce dynamic maps for monitoring the temperature variations in a self-burning coal waste pile, under a GIS open source environment—GIS-ECOAL (freely available). The performance of the application was evaluated with distributed temperature measurements gathered in the S. Pedro da Cova (Portugal) coal waste pile. In order to obtain the temperature data, an optical fiber cable was disposed over the affected area of the pile, with 42 location stakes acting as precisely-located control points for the temperature measurement. A monthly data set from July (15 min of interval) was fed into the application and a video composed by several layouts with temperature measurements was created allowing for recognizing two main areas with higher temperatures. The field observations also allow the identification of these zones; however, the identification of an area with higher temperatures in the top of the studied area was only possible through the visualization of the images created by this application. The generated videos make possible the dynamic and continuous visualization of the combustion process in the monitored area.


Introduction
Despite the benefits to the economic and social sectors of many countries, both coal mining and coal consumption are responsible for impacts that are harmful to the environment and human health [1,2].Those related to coal mining may include changes in landscape and land use, soil erosion, increased noise generation, production of solid wastes, air pollution, surface and groundwater pollution, soils and sediments pollution, land instability, subsidence, coal mine fires and several impacts on local biodiversity [2].Spontaneous combustion of coal has been recognized all over the world [3][4][5] and may occur during coal mining (opencast and underground), storage, waste disposal and transportation.In Portugal, the occurrence of coal related fires was identified in three coal waste piles resulting from mining of the Douro Coalfield [6,7].Some of the concerns associated with coal fires are the loss of natural resources and the impacts on environment and human health.These impacts are related to the generation of large quantities of products that may be released into the atmosphere, which include the emissions of carbon dioxide, sulfur and nitrogen oxides, volatile organic compounds, and particulate matter [3,8].The type and concentration of compounds emitted to atmosphere largely depends on the combustion temperatures.However, sometimes it is difficult to quantify these factors, so a qualitative analysis is therefore combined with a statistical approach in order to correlate multiple data [9].The surface and subsurface coal fire areas can be detected with different techniques, such as borehole temperature measurement, geophysical methods, thermo-compositional analysis, thermographic measurements, airborne remote sensing and space borne remote sensing techniques [10].For instance, Gangopadhyay et al. [11] attempt to identify and extract the temperature anomalies related to coal fires on an experimental basis.Chen et al. [12] proposed an approach of a basis for a large-scale and low-cost topographic survey for sustainable environmental planning and for mitigating the environmental anthropogenic impacts due to mining, using Unmanned Aerial Vehicles (UAVs) and digital terrain analysis based on elevation data.Geographical Information Systems (GIS) have also been used in studies related to coal fires and to assess and monitor the environmental impacts caused by mining activities.Some authors, such as Prakash and Vekerdy [13], Voigt et al. [14], Zhang et al. [15], and Prakash et al. [16], focused specifically on GIS applications related to coal fires and in the development of tools for detection and analysis.Ali et al. [17] proposed an integration considering Wireless Sensor Networks (WSNs) assisted GIS to monitor and control underground mining applications from a surface office, promoting safety and health, operational management and cost-effectiveness.Wu and Liu [18] established a regional underground coal fire risk assessment (UCF-RA) index system, focused on obtaining the probability of coal combustion in spatial extent and burn intensity.In order to calculate the factors and indexes affecting coal fire risk, ArcGIS 9.0 (ESRI, Redlands, Califórnia, EUA) was used, considering weight assignment.Ozdeniz et al. [19] referred to the difficulty in evaluating internal reactions in coal waste piles.In their study, temperature changes were evaluated through a special mechanism, and a 3.5 D temporal temperature distribution model at two levels was created to evaluate the temperature distribution and the temperature changes, and to understand the behavior of temperature in an industrial coal stockpile.A 3.5 D model can be constructed to represent the distribution of the attribute data, in this case with respect to the spatial information in 3D space.The coal stockpile was equipped with thermocouple sensors, which were connected to a digital transformer panel and then to a computer, in order to automatically collect the measurements.The location-based information analysis technique was connected to GIS, which provides a 3.5 D model in a virtual environment, with respect to time and location of points.The dataset was stored in GIS database.The data was filtered using a moving average filter and the temperature changes over time were plotted.
GIS tools combined with geostatistical techniques provide relevant information to understand the spatial variability of a certain parameter, such as forest fire risk [20][21][22][23].Interpolation algorithms and geostatistical methods are usually available in GIS software.Gharechelou et al. [24] map the soil moisture using interpolation methods and geostatistical techniques such as kriging.In addition, interpolation methods proved to be able to map important features.Lee and Park [25] analysed the risk of ground subsidence where abandoned underground mines exist, through relevant factors such as slope, geology, land use, among others.A decision tree approach was used under GIS software.Song et al. [26] also focused on the prediction of mining subsidence and its impact on the environment with quantitative GIS analysis.An approach to identify the key factors of the impact on soil erosion induced by coal mining subsidence, based on the Revised Universal Soil Loss Equation (RUSLE) model and using GIS, was also explored by Lei et al. [27,28].Yenilmez et al. [29] used GIS tools to evaluate the pollution levels at an abandoned coal mine site at Ovacik-Yaprakli (Cankiri, Turkey) with respect to topography and surface runoff, using a GIS software.It was demonstrated that GIS analysis tools are useful to identify areas with a high concentration of pollutants, which helps to prevent the overlook of highly contaminated sites.Hannemann et al. [30] searched for a simple methodology to store deformation data combined with visualization and analysis of available height change information, using an ArcGIS extension developed for a German hard coal mining company RAG Deutsche Steinkhole, Herne, Deutschland.The developed extension accesses the spatial data through ArcSDE for Oracle database, metadata and remote sensing data.Oh and Lee [31] used GIS techniques to analyze the ground subsidence hazard in abandoned coal mines in Samcheok, Korea.A Digital Elevation Model (DEM) was used to calculate the slope, depth of drift, distance from drift, groundwater level, permeability, geology and land use from different databases.Using this data, the relations between the ground subsidence and the seven factors referred to before were quantified through frequency ratio, weights-of-evidence, logistic regression and artificial neural network models.They conclude that the integration model is better than the individual model.
All the studies referred to above were developed in GIS proprietary software, which requires a license to use.Based on the literature, there is not a GIS application to evaluate and monitor the temperature of self-burning coal waste piles.In addition, and considering that the management of the coal waste piles with continuous control and identification of probable evolution scenarios through temperature measurement is crucial to minimize the impacts of the process, an application based on an open source GIS software, allowing for evaluation and monitoring temperature in self-burning coal waste piles that would be very valuable in this field of geosciences.Such application would offer the possibility to modify the code and adapt it to any region with different characteristics, providing an added value for coal related fire control and monitoring.
The aim of this paper was to develop a new application to produce dynamic temperature maps allowing the monitoring of temperature variations in a self-burning coal waste pile, under a GIS open source environment.It was also expected to create an application with a user-friendly graphic interface, where the user can introduce the temperature measurements, and obtain several interpolation maps and a video, which allows for evaluating the temperatures dynamic in a specific area.The application was developed within the QGIS open source software, and it is free to use.

Case Study-ECOAL MGT Project (Ecological Management of Coal Waste Piles in Combustion)
The occurrence of coal related fires in Portugal was identified in three coal waste piles resulting from mining activities in Douro Coalfield (NW of Portugal).The Douro Coalfield, from terrestrial Carboniferous (Ghzeliano, Upper Pennsylvanian (Lower Stephanian C)) [32,33], represents the most important coal-bearing deposit in Portugal, with NW-SE alignment, a variable width (30-250 m) and approximately 53 km in length [34,35].The mining activities in Douro Coalfield started in 1795 and went until 1994.The exploited anthracite A coal type was used for power generation [36].
The referred coal waste piles, namely, S. Pedro da Cova, Lomba, and Midões, were emplaced over the old mine sites and adjacent areas and started burning in 2005, after forest fires that caused their ignition.The S. Pedro da Cova waste pile, which is the biggest and presently still burning, is located very close to the oldest centre of mining activities in S. Pedro da Cova and close to a residential area and social infrastructure facilities.This pile has an elongated form, occupying an area of approximately 28,000 m 2 , with the burning process occurring only in the south slope and moving along it.Figure 1 shows the area that is currently under the self-burning process.
Previous studies about S. Pedro da Cova waste pile dealt with petrographic, mineralogical and geochemical characterization of the waste material, including the potential leaching of inorganic constituents, and with the identification of gases emitted from gas vents [6,[37][38][39][40].The analysis of gases exhibited a wide range of hydrocarbons, including aromatic and aliphatic hydrocarbons [6], which are known to have deleterious environmental and human health impacts.The types and concentrations of gases emitted during the coal wastes combustion depend, among other factors, on the combustion temperature.It was also demonstrated that the burning coal waste piles may have reached combustion temperatures of about 1000 • C or even higher [39,40].In this context, the monitoring of combustion temperature and gaseous emissions in S. Pedro da Cova waste pile is essential to identify the associated hazards and probable evolution scenarios, allowing for the definition of timely corrective measures to minimize the negative impacts caused by this process.

Multi-Point Measurement Sensing System
Owing to the extreme conditions associated with the high temperatures caused by the combustion process, the multi-point measurement sensing system was implemented employing optical fiber sensor technology.Fiber sensors have unique features, among them their capability to measure in harsh environments, high sensitivity, durability, and remote measurement capabilities.In particular, a distributed optical fiber sensor was utilized to monitor the combustion temperature in the coal waste pile.This kind of sensor allows for extracting the values of temperature as a function of position along the length of a sensing fiber cable.This feature is particularly attractive for use in applications where monitoring of temperature is required over a large number of points covering a wide area, since the whole fiber acts as a sensor, replacing a vast number of discrete sensors with a single optical fiber cable.
Distributed optical fiber sensing techniques are based on some kind of light scattering processes taking place inside the optical fiber.In this case, Stimulated Brillouin Scattering (SBS) is at the heart of the employed technique, denominated Brillouin Optical Time Domain Analysis (BOTDA) [41].The architecture of the BOTDA prototype used in S. Pedro da Cova is similar to a conventional BOTDA [42], but incorporating novel techniques to improve the spatial resolution.In this case, a differential pulse-width pair technique [43] was used to obtain a spatial resolution of ~0.5 m, and the signal-tonoise ratio of the measurements was enhanced by using balanced detection [44] and polarization diversity [45].
For the temperature distributed measurement, ~1 km of a special dual-fiber optical cable that can withstand high temperature was extended over the whole area of the waste pile in combustion (BRUsens Temperature sensing cable 5.0 mm High-density polyethylene (HDPE), provided by Brugg Cables).Prior to the installation, the cable was calibrated in the lab so as to have a reference of Brillouin Frequency Shift (BFS) distribution and sensitivity at each position along the fiber cable.Figure 2 shows a calibration of sensitivity at one particular position along the fiber cable in the range of 5-218 °C.As it is visible, the response of the BFS with temperature of the fiber used is essentially linear with typical sensitivity values in the order of 1.1 MHz/°C.

Multi-Point Measurement Sensing System
Owing to the extreme conditions associated with the high temperatures caused by the combustion process, the multi-point measurement sensing system was implemented employing optical fiber sensor technology.Fiber sensors have unique features, among them their capability to measure in harsh environments, high sensitivity, durability, and remote measurement capabilities.In particular, a distributed optical fiber sensor was utilized to monitor the combustion temperature in the coal waste pile.This kind of sensor allows for extracting the values of temperature as a function of position along the length of a sensing fiber cable.This feature is particularly attractive for use in applications where monitoring of temperature is required over a large number of points covering a wide area, since the whole fiber acts as a sensor, replacing a vast number of discrete sensors with a single optical fiber cable.
Distributed optical fiber sensing techniques are based on some kind of light scattering processes taking place inside the optical fiber.In this case, Stimulated Brillouin Scattering (SBS) is at the heart of the employed technique, denominated Brillouin Optical Time Domain Analysis (BOTDA) [41].The architecture of the BOTDA prototype used in S. Pedro da Cova is similar to a conventional BOTDA [42], but incorporating novel techniques to improve the spatial resolution.In this case, a differential pulse-width pair technique [43] was used to obtain a spatial resolution of ~0.5 m, and the signal-to-noise ratio of the measurements was enhanced by using balanced detection [44] and polarization diversity [45].
For the temperature distributed measurement, ~1 km of a special dual-fiber optical cable that can withstand high temperature was extended over the whole area of the waste pile in combustion (BRUsens Temperature sensing cable 5.0 mm High-density polyethylene (HDPE), provided by Brugg Cables).Prior to the installation, the cable was calibrated in the lab so as to have a reference of Brillouin Frequency Shift (BFS) distribution and sensitivity at each position along the fiber cable.Figure 2 shows a calibration of sensitivity at one particular position along the fiber cable in the range of 5-218 • C. As it is visible, the response of the BFS with temperature of the fiber used is essentially linear with typical sensitivity values in the order of 1.1 MHz/ • C. The cable was buried at 10 to 20 cm depth.Two optical fibers were used in the measuring cable, so redundant measurements were obtained in each position.The temperature estimation was therefore achieved by averaging the results obtained from the two independent measurements realized in the same position [46].The BOTDA interrogation unit was placed in a data collection station situated close to the coal waste pile in combustion and a standard optical cable (~500 m) was utilized in order to connect the special optical cable with the BOTDA system according to the schematic representation in Figure 3.A comprehensive set of temperature measurements in the coal waste pile were accomplished over several months by employing this BOTDA sensor, contributing to the identification of the associated risks, the prediction of the combustion dynamics, as well as the establishment of evolution scenarios.
The optical fiber cable was installed with 42 location stakes, in order to support the control points.According to the field conditions, the optical fiber was extended in a straight line between two consecutive points.In the other points, the optical fiber was extended in a straight line, and, at The cable was buried at 10 to 20 cm depth.Two optical fibers were used in the measuring cable, so redundant measurements were obtained in each position.The temperature estimation was therefore achieved by averaging the results obtained from the two independent measurements realized in the same position [46].The BOTDA interrogation unit was placed in a data collection station situated close to the coal waste pile in combustion and a standard optical cable (~500 m) was utilized in order to connect the special optical cable with the BOTDA system according to the schematic representation in Figure 3.The cable was buried at 10 to 20 cm depth.Two optical fibers were used in the measuring cable, so redundant measurements were obtained in each position.The temperature estimation was therefore achieved by averaging the results obtained from the two independent measurements realized in the same position [46].The BOTDA interrogation unit was placed in a data collection station situated close to the coal waste pile in combustion and a standard optical cable (~500 m) was utilized in order to connect the special optical cable with the BOTDA system according to the schematic representation in Figure 3.A comprehensive set of temperature measurements in the coal waste pile were accomplished over several months by employing this BOTDA sensor, contributing to the identification of the associated risks, the prediction of the combustion dynamics, as well as the establishment of evolution scenarios.
The optical fiber cable was installed with 42 location stakes, in order to support the control points.According to the field conditions, the optical fiber was extended in a straight line between two consecutive points.In the other points, the optical fiber was extended in a straight line, and, at A comprehensive set of temperature measurements in the coal waste pile were accomplished over several months by employing this BOTDA sensor, contributing to the identification of the associated risks, the prediction of the combustion dynamics, as well as the establishment of evolution scenarios.
The optical fiber cable was installed with 42 location stakes, in order to support the control points.According to the field conditions, the optical fiber was extended in a straight line between two consecutive points.In the other points, the optical fiber was extended in a straight line, and, at the end of the lines, the transition to the next line was defined by a semi circumference (Figure 3, represented as yellow dots).
In the optical fiber, the temperature values were measured by the BOTDA sensor every 20 cm and spatially referenced with a number and the distance to the origin.The stakes were positioned coinciding with marks printed in the coating of the optical cable to indicate the distance to the origin.A Global Positional System (GPS) survey was performed using a dual frequency Trimble R6 Global Navigation Satellite System (GNSS) receiver, with real time corrections provided by a national network of GPS permanent stations network (RENEP-Rede Nacional de Estações Permanentes GNSS), operated by Direcção Geral do Território [47].The method adopted ensures an accuracy of 2 to 3 cm, in 3D coordinates, which is adequate for the study case.Coordinates of the stakes were calculated in this way for all the stakes, in the national reference system (ETRS89/PTTM06).According to the spatial arrangement of the fiber, the coordinates of 3281 sensor locations were estimated, and stored in the shapefile format, for further analysis in a GIS.A monthly data set from July (with 15 min of interval between measurements) was evaluated and processed.July was selected due to the lack of continuous data in other months.

Open Source GIS Software and Programming Language
The application described in this work was developed under the GIS open source software, QGIS version 2.10 Pisa (Manufacturer, City, US State abbrev.if applicable, Country?).QGIS is a free and open source software, subject to the GNU General Public license.QGIS was chosen due to the authors' previous experience and skills [22,23,48,49].It has a user-friendly graphic interface and several robust algorithms/plugins to view, edit, manipulate and analyze geospatial data.QGIS is developed in C++ programming language, but it is complemented with plugins developed in Python [50].It has the advantage of being integrated with other open source GIS packages, including PostGIS, GRASS, SAGA, GDAL, and Orfeo Toolbox, among others.In the scope of this work, QGIS API, QT API, PyQt API, with QtCore, QtGui, Phonon and GDAL were some of the libraries used [51,52].QGIS methods complemented with Python methods allow for the creation of different multimedia formats, which, in the scope of this work, can be very useful to understand the dynamics of the temperature evolution.
Python was used in the development of this application.It is free and open source, developed under an Open Source Initiative (OSI)-approved license, administrated by the Python Software Distribution.It is friendly to use and learn and users have the support of very useful documentation and mailing lists [50].Several modules from Python were used in this work: re, glob, subprocess, os and sys are some of them.

Graphic Interface
The GIS-ECOAL application is composed of a single button, which can be added to QGIS.It opens a window consisting of eight fields divided in three sections (Figure 4).The first section, called "Input temperature data and study zone information", allows for the input of temperature measurement data, in text format, and boundaries of the study area, in Environmental Systems Research Institute (ESRI) shapefile format.The next section ("Input coordinates information") is to input the location of sampling points, which can be provided in the form of an xy text file, or as a point shapefile.The third section is to specify the path for the raster output data, in TIFF format, which will be used to compose a video.The video will be displayed in a window, initially presented in black.A plot can be generated based on the temperature measurements for a specific point (defined as shapefile) through the "Plot" button.Figures 4 and 5 present the user case diagram and the application graphic interface, respectively.The application was developed considering widgets such as combo boxes, buttons and video widgets.The input and output folders were connected through QDir class of Qt API, which allows for having access and manipulating directory structures and their contents [53].In order to develop the GIS-ECOAL application, two scripts were created, containing the data manipulation and the graphic interface definition. Figure 7 presents the application class diagram.The application starts to verify the text file extension, and, if there are strings, their content.In addition, the temperature data is given with float format.It was defined that the relevant position values, for this specific work, range between 240 m and 935 m approximately, taking into account that the position data is given starting in zero meters.In order to restrict the values to these intervals, an if condition was created and only the values in that specified interval were written.Then, the values were grouped in ten elements and the median value of these ten elements was calculated, using methods as sorted, int, len and round.Only one point is retained in every 10, so that data becomes filtered of some remaining noise and the point density along the fibre becomes similar to the density across.This is important in order to avoid biasing in the interpolation methods described below.
As the application gives two options in order to introduce the point coordinates of the study zone as input, an if statement was created.In the first case, if a text file with the coordinates disposed in two columns, x and y, is introduced, this input file (Input TXT coordinate points) is used to the The application was developed considering widgets such as combo boxes, buttons and video widgets.The input and output folders were connected through QDir class of Qt API, which allows for having access and manipulating directory structures and their contents [53].In order to develop the GIS-ECOAL application, two scripts were created, containing the data manipulation and the graphic interface definition.The application was developed considering widgets such as combo boxes, buttons and video widgets.The input and output folders were connected through QDir class of Qt API, which allows for having access and manipulating directory structures and their contents [53].In order to develop the GIS-ECOAL application, two scripts were created, containing the data manipulation and the graphic interface definition. Figure 7 presents the application class diagram.The application starts to verify the text file extension, and, if there are strings, their content.In addition, the temperature data is given with float format.It was defined that the relevant position values, for this specific work, range between 240 m and 935 m approximately, taking into account that the position data is given starting in zero meters.In order to restrict the values to these intervals, an if condition was created and only the values in that specified interval were written.Then, the values were grouped in ten elements and the median value of these ten elements was calculated, using methods as sorted, int, len and round.Only one point is retained in every 10, so that data becomes filtered of some remaining noise and the point density along the fibre becomes similar to the density across.This is important in order to avoid biasing in the interpolation methods described below.
As the application gives two options in order to introduce the point coordinates of the study zone as input, an if statement was created.In the first case, if a text file with the coordinates disposed in two columns, x and y, is introduced, this input file (Input TXT coordinate points) is used to the The application starts to verify the text file extension, and, if there are strings, their content.In addition, the temperature data is given with float format.It was defined that the relevant position values, for this specific work, range between 240 m and 935 m approximately, taking into account that the position data is given starting in zero meters.In order to restrict the values to these intervals, an if condition was created and only the values in that specified interval were written.Then, the values were grouped in ten elements and the median value of these ten elements was calculated, using methods as sorted, int, len and round.Only one point is retained in every 10, so that data becomes filtered of some remaining noise and the point density along the fibre becomes similar to the density across.This is important in order to avoid biasing in the interpolation methods described below.
As the application gives two options in order to introduce the point coordinates of the study zone as input, an if statement was created.In the first case, if a text file with the coordinates disposed in two columns, x and y, is introduced, this input file (Input TXT coordinate points) is used to the next step.In the second case, if a shapefile with points is introduced (Input Sahpefile (SHP) points), the application reads the vector file, iterates over the geometry features and extracts the coordinates of each point, through asPoint() function belonging to QgsGeometry [53].In the next step, the coordinates, x and y, are added along with the median temperature values in a text file composed by three columns: median temperature values; x-coordinate; and y-coordinate.The text files were converted to shapefile using the writeAsVectorFormat() function from QgsVectorFileWriter class.With the points created and the temperature values in the attribute table, the ordinary kriging algorithm from SAGA was employed in order to interpolate the points and generate a surface.This algorithm allows for generating a grid interpolation from irregular sample points [54].In order to restrict the surfaces extension to the study area, the clip grid with polygon algorithm from SAGA was used.In order to obtain a graphical representation of the range of colour temperatures in degrees an adequate colour ramp was defined for each surface, using QgsColorRampShader class, a ramp shader which colours a raster based on a list of value ranges in a ramp [53].A colour ramp with green colour for low values of temperature and red colour for high values of temperature was assigned to each surface through QgsColor.The range of values used was the adopted to all the images.Several functions of QGIS API were applied to perform this task, setColorRampItemList, setRasterShaderFunction, QgsSingleBandPseudoColorRenderer class, among others.In order to create a composed layout with the elements, title, graphic scale, legend and north arrow, the QgsComposition for map composer was used.The application reads the coordinate system of the surfaces in map canvas, sets an extent, and loads a predefined template.In this template, the application adds the map (temperature surface), the legend with the colour ramp defined, the title which is assigned according to the name of the file, the north arrow through an image, and the graphic scale bar according to the coordinate system of the surfaces.
Finally, the composed maps were exported in a tif format with the name Temperature (Celsius) and a sequential number and saved in the output folder.
Ffmpeg, which is a single command line program for multimedia format convertion, was used to compose the video from the sequence of individual images.Ffmpeg calls several libraries such as libx264, which is the codec used in this work [55].The command line ffmpeg was installed in the system, and it was called through subprocess module from Python.This module allows for spawning new processes, connect to their input, output or error pipes and obtain their return codes [50].The application creates the video, in webm format, which can be displayed in any video player programme.However, the application provides the possibility of video visualization directly in the application graphic interface.The button Play is connected to the function PlayClicked, which reads the images directory, verifies if the video exists and, if so, loads and plays it in the application.To perform this task, the MediaSource class from Phonon Class Reference belonging to PyQt4 library was used [51].The Phonon namespace is composed for classes and functions for multimedia applications [56].
The user has the possibility to evaluate the temperature dynamic of a specific point in the field, through a graphic representation of the temperature values as function of time and a table with a statistical analysis for each day.In order to generate this plot, an input field (Input point shapefile) and a button (Plot) were created.The first field defines the shapefile with the specific point and the button allows for generating the graphic with the dynamic of the temperature.The maximum and minimum temperature values are presented as text next to the plot.To perform this task, a function, plot_create, was developed and matplotlib library was imported.Matplotlib provides several plotting capabilities with useful modules and functions to handle dates as datetime, date2num, num2date, among others.The figure module from matplotlib was also used in order to define the plot elements.The final plot can be exported to PDF format or image format [57].The large amount of data can be very difficult to represent in the x-axis, so the temperature median value for a day was calculated using Python dictionaries.A txt file is also automatically created through Python libraries.The statistical variables are calculated through numpy library functions.

Experimental Validation
Spatial continuous data is usually required for environmental management and monitoring studies.Temperature data and other climate variables are normally recorded through disperse weather stations and collected from points sources, so the data spatial distribution requires a spatial interpolation in order to provide climate data in other sites [43,58].Several spatial interpolation methods are widely used such as Inverse Distance Weighting (IDW), Delaunay Triangulation and Geostatistics The latter analyses the spatial correlation between the data, using the climate data recorded in the stations or combining this information with geographic or topographic information [58].
Ordinary kriging was used in this work, considering only the data recorded in the stations.Ordinary kriging assumes that the data are stationary, i.e., they contain no significant trends (or drift) over the spatial range [59].Kriging can be described as follows: let x 1 , x 2 , x 3 , . . ., x n be points of observation with observed values Z(x 1 ), Z(x 2 ), Z(x 3 ), . . ., Z(x n ), and x 0 an unsampled point with unknown value Z(x 0 ).The estimation Z*(x 0 ) of that value is obtained using a linear combination of the n known values [60]: Coefficients λ i are determined in such way that the estimation Z* (x 0 ) is unbiased (Equation ( 2)): and that the estimation variance given by Equation ( 3) is minimal.
The conditions given by Equations ( 2) and (3) will be respected if coefficients λ i are the solution of the equation given in Equation ( 4), where where γ is the theoretical semi-variogram and µ is a Lagrange multiplier.The variogram is a function describing the dependence of the variance with respect to distance.In order to apply kriging, a simplified model should be fitted to a variogram.Usually, there are three parameters that describe a variogram with enough detail: nugget, sill and range.More details about semi-variograms, and these parameters can be founded in Cressie [61].Two other interpolation methods-IDW and Triangulation-were also applied, in order to validate the results obtained by kriging.When compared to kriging, IDW has the disadvantage of doing an average of a certain number of neighbour points, without any spatial correlation analysis.The Triangulation method uses the three neighbour points forming a triangle, and tends to also model any noise that the data may have.It also has the disadvantage that it only interpolates inside a convex hull defined by the data points.
The kriging interpolation was performed through ordinarykriging algorithm from SAGA.The IDW interpolation method was performed through v.idw algorithm from GRASS [62].The differences between the measured temperature value and the interpolated ones (residuals) were calculated in order to find the more accurate interpolation method.

Results and Discussion
In order to assess and compare the results of the methods, a residual analysis on the data points was made.A residual is the difference between the point value and the result of the value estimated by the grid, using bilinear interpolation.Then, statistics of the residuals, such as average, standard deviation, Root Mean Square (RMS), minimum and maximum, were calculated.Since, in general, the average of the residuals is close to zero, the standard deviation and the RMS are very similar.The RMS values were also expressed as relative errors, in proportion of the amplitude interval of the temperatures, which was around 78 • C.
The two methods were applied, using the same sample points of optical fiber, in order to generate grids with 1 meter spacing, which is denser than the sampling.
Figure 8 shows, for the three methods, the grid generated and the residuals on the data points, with colours associated to the values.All maps were obtained with QGIS.
Table 1 presents the statistical analysis for each interpolation method.
ISPRS Int.J. Geo-Inf.2017, 6, 87 11 of 19 The two methods were applied, using the same sample points of optical fiber, in order to generate grids with 1 meter spacing, which is denser than the sampling.
Figure 8 shows, for the three methods, the grid generated and the residuals on the data points, with colours associated to the values.All maps were obtained with QGIS.
Table 1 presents the statistical analysis for each interpolation method.From the results obtained, the residuals between −2 • C and 2 • C were recorded in 59% of the samples in the case of kriging, 83% for Triangulation and 48% for the IDW method.Twenty-one percent of the samples presented differences higher than 2 • C through kriging interpolation.Through the IDW interpolator, the percentage of samples with residuals higher than 2 • C was 20% and 10% through Triangulation.With respect to residuals lower than −2 • C, kriging interpolation presents 20% of samples, 32% with IDW and 7% with Triangulation.The statistical results obtained were very similar; however, IDW is a more simplistic method.
Regarding the statistics, the Triangulation method presented the best accuracy, with the smaller RMS.However, this method is restricted to the convex hull defined by the sampling point, so only 305 points could be evaluated.IDW and kriging presented similar statistical results and both can be applied to the whole area.Nonetheless, kriging provides an optimal interpolation estimate for a given coordinate location, as well as a variance estimate for the interpolation value, while IDW estimation is only based on weights given by distance from the interpolation location.Therefore, based on this, the kriging method was chosen to interpolate the data.
In order to test the application described, the S. Pedro da Cova waste pile was evaluated for the study of the temperature dynamics in a certain period of time.To evaluate the application performance, the processing time was recorded.The data set was composed by 782 files, which took approximately three hours to process and to create 782 images in an Intel ® Celeron ® CPU N2820 @ 2.13 GHz (Porto, Portugal), with 4.00 GB of RAM, and a Windows 8.1 operating system (Porto, Portugal).The video was created in approximately one hour.Figure 9 presents two images created at different times, 05h05 and 18h05, on the same day (11 July), where graphical elements, such as north arrow, legend and graphic scale were added.The user can play the resulting video by pressing the Play button.The temperature images presented in the final layouts were obtained by an ordinary kriging interpolation.From the results obtained, the residuals between −2 °C and 2 °C were recorded in 59% of the samples in the case of kriging, 83% for Triangulation and 48% for the IDW method.Twenty-one percent of the samples presented differences higher than 2 °C through kriging interpolation.Through the IDW interpolator, the percentage of samples with residuals higher than 2 °C was 20% and 10% through Triangulation.With respect to residuals lower than −2 °C, kriging interpolation presents 20% of samples, 32% with IDW and 7% with Triangulation.The statistical results obtained were very similar; however, IDW is a more simplistic method.
Regarding the statistics, the Triangulation method presented the best accuracy, with the smaller RMS.However, this method is restricted to the convex hull defined by the sampling point, so only 305 points could be evaluated.IDW and kriging presented similar statistical results and both can be applied to the whole area.Nonetheless, kriging provides an optimal interpolation estimate for a given coordinate location, as well as a variance estimate for the interpolation value, while IDW estimation is only based on weights given by distance from the interpolation location.Therefore, based on this, the kriging method was chosen to interpolate the data.
In order to test the application described, the S. Pedro da Cova waste pile was evaluated for the study of the temperature dynamics in a certain period of time.To evaluate the application performance, the processing time was recorded.The data set was composed by 782 files, which took approximately three hours to process and to create 782 images in an Intel ® Celeron ® CPU N2820 @ 2.13 GHz (Porto, Portugal), with 4.00 GB of RAM, and a Windows 8.1 operating system (Porto, Portugal).The video was created in approximately one hour.Figure 9 presents two images created at different times, 05h05 and 18h05, on the same day (11 July), where graphical elements, such as north arrow, legend and graphic scale were added.The user can play the resulting video by pressing the Play button.The temperature images presented in the final layouts were obtained by an ordinary kriging interpolation.Images relative to selected days and specific times can be retrieved.The analysis of individual images gives valuable information about the variation of the temperature in different moments of a day (i.e., as shown in Figure 9) and also about the distribution and dispersion of the combustion activity and intensity in the fire area, since the combustion process in a coal waste pile is spatially very heterogeneous.
These images also allow the identification of the temperatures associated to the combustion process and the main focus and/or new focus of the combustion activity.Taking Figure 9 as an example, one main area of more intense combustion can be recognized, which is the large red part at the bottom of the monitored area.Figure 9 also shows that the temperature of the waste pile material is generally higher at 18:05 of 11th of July when compared with the temperatures measured at 05h05, which is related to the influence of the atmospheric conditions, which will be discussed later.Field observations and thermographic analysis of the studied area had allowed the identification of this main focus of the combustion activity; however, the identification of significantly higher temperatures in the left and upper corner of the waste pile (more yellowish areas) indicate that the combustion activity may be increasing in this area.The identification of this new focus of combustion was only possible through the visualization of the images created by the GIS-ECOAL application.The above points out the significance of the developed application and the need for continuous temperature monitoring with sensing systems.
The integration of images relative to a given period of time (a day, a week, a month, etc.) provides information about the variation of the temperature due to atmospheric conditions and seasonal variations.The generated videos make possible the dynamic and continuous visualization of the combustion process in the monitored area.The data provided by the GIS-ECOAL application allows: (i) the identification of new combustion areas; (ii) the monitoring of combustion focus; and, (iii) the prediction of evolution scenarios that may contribute to the establishment of mitigation measures and appropriate management practices.
In order to test the functionality to create a plot for a specific point and represent the temperature evolution in time, three points were chosen according to the video analysis.The location of the selected points for detailed analysis can be observed in Figure 9: point 1 in the main focus of combustion (also recognized in the field); point 2 in a potential new focus of combustion; and point 3 in an area without combustion activity.For each point, a shapefile was generated and introduced as an input in the application.The application runs through the images created and assigns the median temperature value in each day and time of the point.The values were saved, and, finally, a plot as a function of time was created.The median value was chosen due to the large amount of data.In addition, a text file is automatically created with the statistical variables related to each day, such as average, median value, maximum and minimum temperature, standard deviation, coefficient of variation and amplitude.This statistical analysis allows the data evaluation for each day.Figure 10 represents the resulting plots for points 1, 2, and 3, which illustrate the temperature dynamics along July 2015.The x-axis shows the evaluated days, and the temperature measurements are presented on the y-axis, as the median for each day.The plots are generated for visualization within QGIS, but the user can save them as an image or PDF, for evaluation at any other time.
The analysis of Figure 10 demonstrates that temperatures at point 1, which is located in an area with intense combustion activity, are higher than 90 • C while, at point 2, the temperatures vary between 50  The atmospheric conditions can strongly influence the temperature of the solid materials that constitute the coal waste pile, as well as the combustion temperature.Based on the data from the Pedras Rubras meteorological station (41°9′0″ N, 8°37′12″ W), the variation of temperature and humidity at points 1, 2 and 3, in the hottest and coldest days of July (considering the maximum temperatures) (Table 2), were evaluated.The differences of temperature and humidity between the coldest and hottest days are not very significant, as July is a summer month in the northern hemisphere (Table 2).Figure 11 presents the graphical representation of temperature variation at points 1, 2, and 3 and the representation of atmospheric temperature of the coldest and hottest days of July.It can be observed that the maximum atmospheric temperature of the analysed days and the maximum combustion temperature of coal waste pile at the three points did not occur at the same time of the day.The maximum temperature at points 1, 2, and 3 was reached about three to five hours later than the maximum atmospheric temperature.The same is observed for the minimum temperature of those days when compared with the atmospheric temperature.The displacement of the temperature of coal combustion of the waste pile and the atmosphere is attributed to the progressive heat transfer between the atmosphere and solid materials.In addition, at points 1 and 3, the variation of the temperature along the days is more irregular than the temperature at point 2. In the case of point 1, this fact could be attributed to the occurrence and duration of the combustion process.The combustion at point 1 was identified and has been monitored since 2013 through thermographic images obtained once a month and normally by midday [39].The atmospheric conditions can strongly influence the temperature of the solid materials that constitute the coal waste pile, as well as the combustion temperature.Based on the data from the Pedras Rubras meteorological station (41 • 9 0" N, 8 • 37 12" W), the variation of temperature and humidity at points 1, 2 and 3, in the hottest and coldest days of July (considering the maximum temperatures) (Table 2), were evaluated.The differences of temperature and humidity between the coldest and hottest days are not very significant, as July is a summer month in the northern hemisphere (Table 2).Figure 11 presents the graphical representation of temperature variation at points 1, 2, and 3 and the representation of atmospheric temperature of the coldest and hottest days of July.It can be observed that the maximum atmospheric temperature of the analysed days and the maximum combustion temperature of coal waste pile at the three points did not occur at the same time of the day.The maximum temperature at points 1, 2, and 3 was reached about three to five hours later than the maximum atmospheric temperature.The same is observed for the minimum temperature of those days when compared with the atmospheric temperature.The displacement of the temperature of coal combustion of the waste pile and the atmosphere is attributed to the progressive heat transfer between the atmosphere and solid materials.In addition, at points 1 and 3, the variation of the temperature along the days is more irregular than the temperature at point 2. In the case of point 1, this fact could be attributed to the occurrence and duration of the combustion process.The combustion at point 1 was identified and has been monitored since 2013 through thermographic images obtained once a month and normally by midday [39].

Conclusions
The continuous monitoring of waste piles in self-combustion is relevant in order to minimize environmental impacts.The application developed in this work allows for the monitoring of combustion temperature in a self-burning coal waste through the creation of images and videos representing the spatial and temporal variations of temperature, as well as doing graphical analysis in time for specific points.To test the application, the case study of the S. Pedro da Cova coal waste pile in Portugal was used.The processing time was recorded in order to evaluate the application performance, which proved to be acceptable.The results allow for the interpretation and analysis of the dynamics of the temperature and combustion process at specific periods of time.The graphic representation of temperature for certain points was performed with the median temperature value for selected days due to the time distribution of data.The resulting plots can be very useful to analyse and evaluate the behaviour of the combustion process at specific points.The developed application is easy to use and manipulate, even for a user not familiar with the QGIS environment.The information obtained from the GIS-ECOAL application can be particularly useful to predict evolution

Conclusions
The continuous monitoring of waste piles in self-combustion is relevant in order to minimize environmental impacts.The application developed in this work allows for the monitoring of combustion temperature in a self-burning coal waste through the creation of images and videos representing the spatial and temporal variations of temperature, as well as doing graphical analysis in time for specific points.To test the application, the case study of the S. Pedro da Cova coal waste pile in Portugal was used.The processing time was recorded in order to evaluate the application performance, which proved to be acceptable.The results allow for the interpretation and analysis of the dynamics of the temperature and combustion process at specific periods of time.The graphic representation of temperature for certain points was performed with the median temperature value for selected days due to the time distribution of data.The resulting plots can be very useful to analyse and evaluate the behaviour of the combustion process at specific points.The developed application is easy to use and manipulate, even for a user not familiar with the QGIS environment.The information obtained from the GIS-ECOAL application can be particularly useful to predict evolution scenarios, which, in turn, may contribute to the identification of new combustion focus, mitigation measures and appropriate management practices.
An added value of the application is that it can be adapted and applied to other situations related to self-heating and self-combustion of coal in seams, during mining, handling and transportation and storage piles.The self-heating of landfills represents another environmental problem where the application can be adapted and used.
The desktop GIS application, with all the source code that can be modified or adapted, is available at [63] so that the reader can test the GIS-ECOAL application.

19 Figure 1 .
Figure 1.General view of the active self-burning area in S. Pedro da Cova waste pile (base map source: Bing Aerial, Imagery of NASA Earthstar Geographics SIO© 2016).

Figure 1 .
Figure 1.General view of the active self-burning area in S. Pedro da Cova waste pile (base map source: Bing Aerial, Imagery of NASA Earthstar Geographics SIO© 2016).

Figure 2 .
Figure 2. BFS dependence on temperature for the special optical cable used in S. Pedro da Cova.

Figure 3 .
Figure 3. General view of the active combustion area in S. Pedro da Cova waste pile and schematic representation of the optical fiber system installation in the field (base map source: Google Earth ©2016 Google Image Landsat Data SIO, National Oceanic and Atmospheric Administration (NOAA), U.S. Navy, National Geospatial-Intelligence Agency (NGA), General Bathymetric Chart of the Oceans (GEBCO)).

Figure 2 .
Figure 2. BFS dependence on temperature for the special optical cable used in S. Pedro da Cova.

Figure 2 .
Figure 2. BFS dependence on temperature for the special optical cable used in S. Pedro da Cova.

Figure 3 .
Figure 3. General view of the active combustion area in S. Pedro da Cova waste pile and schematic representation of the optical fiber system installation in the field (base map source: Google Earth ©2016 Google Image Landsat Data SIO, National Oceanic and Atmospheric Administration (NOAA), U.S. Navy, National Geospatial-Intelligence Agency (NGA), General Bathymetric Chart of the Oceans (GEBCO)).

Figure 3 .
Figure 3. General view of the active combustion area in S. Pedro da Cova waste pile and schematic representation of the optical fiber system installation in the field (base map source: Google Earth ©2016 Google Image Landsat Data SIO, National Oceanic and Atmospheric Administration (NOAA), U.S. Navy, National Geospatial-Intelligence Agency (NGA), General Bathymetric Chart of the Oceans (GEBCO)).

Figure 4 .
Figure 4. User case diagram of a GIS-ECOAL application.

3. 3
.2. Application Architecture GIS-ECOAL imports several libraries and Application Programming Interfaces (APIs) in order to create the graphic interface and connections to the different files.

Figure 6
presents the application package diagram.

Figure 4 .
Figure 4. User case diagram of a GIS-ECOAL application.

Figure 4 .
Figure 4. User case diagram of a GIS-ECOAL application.

3. 3
.2. Application Architecture GIS-ECOAL imports several libraries and Application Programming Interfaces (APIs) in order to create the graphic interface and connections to the different files.

Figure 6
presents the application package diagram.

3. 3
.2. Application Architecture GIS-ECOAL imports several libraries and Application Programming Interfaces (APIs) in order to create the graphic interface and connections to the different files.
Figure 6 presents the application package diagram.

Figure 9 .
Figure 9. Temperature image map compositions for 11th July at 05h05 and 18h05.Figure 9. Temperature image map compositions for 11th July at 05h05 and 18h05.

Figure 9 .
Figure 9. Temperature image map compositions for 11th July at 05h05 and 18h05.Figure 9. Temperature image map compositions for 11th July at 05h05 and 18h05.
• C and 60 • C and, at point 3, are generally between 40 • C and 50 • C.

Figure 10 .
Figure 10.Graphical representation of temperatures for points 1, 2 and 3 as a function of time.

Figure 10 .
Figure 10.Graphical representation of temperatures for points 1, 2 and 3 as a function of time.

Figure 11 .
Figure 11.Graphical representations of temperatures at points 1, 2, 3 and atmospheric temperature on the coldest and hottest days.

Figure 11 .
Figure 11.Graphical representations of temperatures at points 1, 2, 3 and atmospheric temperature on the coldest and hottest days.

Table 1 .
Residuals analysis for the interpolation methods used.

Table 1 .
Residuals analysis for the interpolation methods used.Inverse Distance Weighting.2RootMeanSquare.3PercentualRMS represents the RMS relative to the amplitude of the temperatures. 1

Table 2 .
Record of meteorological conditions.

Table 2 .
Record of meteorological conditions.