Generating annual estimates of forest fire disturbance in Canada: the National Burned Area Composite

Additional keywords: boreal forest, fire perimeters, fire refugia, Landsat, NBR, post-fire, unburned islands.

Timely annual statistics on areas affected by fires are critical to understanding, monitoring and modelling carbon cycles and climate change feedbacks and to supporting sustainable forest management (Bernier et al. 2012). There has been considerable interest in methods of mapping and quantifying area burned for monitoring the occurrence, frequency and extent of fires across forested ecosystems (French et al. 2008; Kasischke et al. 2011; Chu and Guo 2014). Operational methods of mapping burned areas include ground-based or airborne global positioning system (GPS) surveys and interpretation of post-fire aerial photography or satellite imagery (Parisien et al. 2006; Kolden and Weisberg 2007; Zell and Kafka 2012). Across and within years, however, there are data quality issues from variation in how operational agencies in Canada map burned areas (Parisien et al. 2006; Hanes et al. 2019; Table S1 available as Supplementary material). This mapping variation results in two problems: the degree to which an accurate estimate of area burned is determined, and its uncertainty due to the data source and mapping method employed at specific fires.

Unburned areas and those minimally affected by fire form unburned islands (i.e. refugia) within a burn polygon perimeter, in combination with water bodies, contribute to overestimation of area burned (Meddens et al. 2018). Better mapping of variability within fires has become increasingly important (Meddens et al. 2016), not only for improving the accuracy of burned area estimates, but also for understanding the effect of fire on ecosystem processes (Lentile et al. 2006; Burton et al. 2008; French et al. 2008; Robinson et al. 2013; Meigs and Krawchuk 2018). In Canada, most provincial and territorial governments, through designated fire management agencies, generate annual burned area products as a set of vector polygons at the end of each fire season to meet their mandate for yearly reporting and management of fire activity (Parisien et al. 2006). This information has been compiled into a national dataset, originally known as the Canadian Large Fire Database representing large fires (≥200 ha) (Stocks et al. 2003). The Canadian National Fire Database (CNFDB) consists of forest fire data that include fire perimeters (polygon data) and fire locations (point data) as provided by Canadian fire management agencies that encompass provinces, territories, and Parks Canada (Natural Resources Canada 2019a). Both datasets are maintained by the Canadian Wildland Fire Information System (CWFIS) of Natural Resources Canada (Lee et al. 2002; Natural Resources Canada 2019a). The concern with these datasets is that much of the burned area reported includes fire refugia and water.

Uncertainty in the estimates of area burned arises because agency data have been derived from a variety of methods and data sources (Table S1, Supplementary material), in combination with inconsistent mapping of unburned areas within fire perimeters (Amiro et al. 2001; Stocks et al. 2003; Parisien et al. 2006; Kolden et al. 2012). Remotely sensed data can be used to reduce this uncertainty, but not with definitive clarity because there are many burned-area mapping products that have been generated across Canada (Fraser et al. 2004; Goetz et al. 2006; Guindon et al. 2014, 2017; White et al. 2017; Guindon et al. 2018). These diverse products apply different detection algorithms, vary in their spatial resolution and spatiotemporal coverage, and differ from the vector polygons generated by the agencies. Reducing uncertainty in the estimates of area burned is needed for studies of trends in wildfire activity (Short 2015). In addition, carbon flux estimates can vary considerably depending on the spatial resolution of the satellite data used (Mascorro et al. 2015) and the assumed confidence of the product used in modelling (Metsaranta et al. 2017).

In Canada, the National Burned Area Composite (NBAC) was first put into operation in 2004 to undertake burned area mapping as an annual exercise performed the year after the fires occur, using the best data source available for each fire event in that season (Natural Resources Canada 2019b, 2019e). The system fulfils the timing requirements of the National Forest Carbon Monitoring, Accounting and Reporting System (NFCMARS), which Canada uses to meet international reporting commitments on its emission and removal of greenhouse gases (GHGs) (Kurz and Apps 2006; Stinson et al. 2011; Metsaranta et al. 2017). NBAC is a fire perimeter map comprising a vector representation of the outer edge of the burned area, with unburned vegetated island features (refugia) and water removed, from which to determine the actual area burned and its geographic location. Representation of the burned area as a vector polygon results in a product that is comparable with those generated by fire management agencies in Canada, while addressing the problem of overestimating area burned.

The purpose of the present paper is to describe the components of the NBAC and the process involved in generating the NBAC as an annual, national polygon-based burned area product for Canada. The paper also presents results from the NBAC for the years 2004 to 2016 to illustrate its outputs, describe and demonstrate the approach used to remove fire refugia, and to provide a first-ever measure of uncertainty, expressed as confidence intervals for the area burned annually in Canada.

The components for generating the NBAC include the data sources, mapping methods and a production process (Fig. 1). The NBAC provides burned area polygons for Canada’s forests, describing where a fire has occurred and how much area has burned, the cause of the fire, the fire’s start and end dates, the data source and the mapping method. It was developed jointly by the Canada Centre for Mapping and Earth Observation (CCMEO) and the Canadian Forest Service of Natural Resources Canada. As a system, the NBAC integrates a suite of post-fire products from agency and federal programs to produce annual composite maps of burned area polygons, based on the best available data at the time of generation. Data are derived from three specific burned area products. The first of these is the CNFDB, which is the source of agency fire polygon data. In areas where these fire products are of coarse spatial detail, maps from finer-spatial-resolution Landsat imagery processed through the Multi-Acquisition Fire Mapping System (MAFiMS; Fig. 2) form the second burned area product described in this paper. The third burned area product, which is based on SPOT VEGETATION/Proba-V coarse-spatial-resolution imagery and the Hotspot and Normalized Difference Vegetation Index (NDVI) Differencing Synergy (HANDS) algorithm (Fraser et al. 2000), provides burned area polygons for fires not reported by the agencies or mapped using MAFiMS.

Fig. 1.  The National Burned Area Composite (NBAC) system procedures to produce a national burned area product with attribution from fire hotspots and Agency, Multi-Acquisition Fire Mapping System (MAFiMS), and Hotspot and Normalized Differenced Vegetation Index Differencing Synergy (HANDS) burned area polygons. NFIREID, National Fire Identifier.

Fig. 2.  Multi-Acquisition Fire Mapping System procedures to map fire events into burned area polygons with attribution using fire hotspots, pre- and post-burn Landsat imagery and ancillary coverage features. For reference, the flowchart is structured as follows: (a) image acquisition and pre-processing; (b) detection and mapping of burn events; and (c) vector polygon extraction. dNBR, differenced Normalized Burn Ratio; NDVI, Normalized Differenced Vegetation Index.

The CNFDB fire polygons are mapped annually by the provincial, territorial and Parks Canada fire management agencies (Table S1: Agency mapping of burned area). The maps contain polygons with descriptor attributes for each fire, such as fire size (hectares), the fire’s start and end dates, the fire’s cause, data source and mapping method. The total burned areas mapped by the various fire management agencies for a given year (Table S1) were used to determine the proportions of burned area by source and mapping method nationally (Table 1). Across and within years, however, there is variation in how operational agencies in Canada map burned areas (Table 1). This variation is attributed to no single post-fire mapping method or data source being common to all fire management agencies across Canada, and more than one approach often being used within a given year (Table S1a to S1l). Data sources and mapping methods influence the quality of the polygon delineation, which creates uncertainties in burned area estimates (Parisien et al. 2006; Kolden and Weisberg 2007). For these reasons, the NBAC considers agency resources as potential sources of data that may or may not be selected in the final compilation for each year.

Table 1.  Proportions of burned area mapped by various data sources and mapping methods across Canada, 2004 to 2016
The area burned estimates were derived from the Canadian National Fire Database polygon or point data available from the Canadian Wildland Fire Information System (http://cwfis.cfs.nrcan.gc.ca/). Values less than 0.1 are presented as 0.1. GPS, global positioning system; MODIS, Moderate Resolution Imaging Spectrophotometer

The MAFiMS was developed to facilitate processing of Landsat imagery to generate more precise delineations of burned area around the fire scar while also excluding unburned forest islands (i.e. fire refugia) and water features. The primary components of this system include image acquisition and pre-processing, detection and mapping of burn events, and extraction of vector polygons (Fig. 2). The image-gathering process begins with identification of the Landsat image Worldwide Reference System path and row where fires have occurred (Fig. 2a). The date of post-burn image acquisition is often during the latter part of the fire season. If cloud-free scenes are unavailable, post-burn images are selected in early spring of the following year, after snowmelt. Pre-burn images are collected 1 year before the fire event, during the same acquisition period as the post-burn image, to ensure similar solar illumination and vegetation phenology (Key 2005). If data from Landsat-7 Enhanced Thematic Mapper with Scan Line Corrector turned off are selected (Fig. S1a Supplementary material), then a Geospatial Data Abstraction Library gap-fill algorithm (Warmerdam 2008) is applied to interpolate missing image values to render the image suitable for burned area mapping (Fig. S1b).

After the Landsat images are downloaded, a top-of-atmosphere reflectance correction (Chander et al. 2009) is applied to both pre- and post-burn imagery, followed by a topographic normalisation procedure (Civco 1989; Law and Nichol 2004) that uses Canadian digital elevation data at a scale of 1 : 250 000 to reduce the topographic effects in high-relief areas (Natural Resources Canada 2000). Next, cloud and water exclusion masks are automatically created from statistical interpretation and thresholding of the red and near-infrared image bands, and agricultural activities are removed by selection of agriculture areas from national land cover classifications (Fisette et al. 2006; Huffman et al. 2006) (Fig. 2b).

The MAFiMS uses an iterative process to map burned area across multiple fire events occurring within the common footprint of the Landsat image pair. Phenological change within unburned vegetation is initially quantified to ensure that it is excluded from the magnitude of spectral changes of areas considered to have been burned. The upper 25th percentile of the distributions of NDVI values in both the pre- and post-burn images defines a region that is likely unburned for purposes of quantifying phenological change. Within this region, the mean and standard deviation of dNBR are computed to define a ‘pseudo-mask’ of unburned vegetation (Fig. 2b). A high-change layer (HCL) threshold is then computed at three standard deviations above the mean dNBR (differenced Normalized Burned Ratio) (French et al. 2008), where dNBR pixel values larger than the HCL threshold become the ‘HCL mask’, an interim product used to identify broad areas within which actual burned pixels are likely to be found. The HCL mask is then converted to polygons, buffered at 1500 m, and filtered with spatially coincident MODIS (Moderate Resolution Imaging Spectroradiometer) hotspots (USDA Forest Service 2020) to create ‘HCL processing units’ (Fig. 2b) that fully encompass individual fire events. The concept of an HCL processing unit reflects the unique conditions of vegetation type, conditions of burn and time of year associated with a particular fire. Depending on when a given fire starts and ends within the mapping year and whether the pre-fire image is collected more than 1 year before the date of the post-fire image, multiple burns and change events (e.g. harvesting, mortality) may be captured within the HCL mask. The sequence of steps from dNBR to verification of several HCL processing units is illustrated in Fig. 3. Multiple burned areas are visible in the pre- and post-burn Landsat images, of which only a portion burned in 2013, the year of interest for this analysis (Fig. 3a and 3b). Several fires are visible because of the magnitude of their spectral changes within the dNBR image (Fig. 3c), and these were subsequently identified in the HCL mask (Fig. 3d). Hotspots are used in the attributes for the start and end dates of the burn, and intersection of hotspots with the HCL mask (Fig. 3e) ensures that only those burns associated with the 2013 hotspots will be used as processing units to map individual fire events (Fig. 3f).

Fig. 3.  An example of multiple burns on a portion of a Landsat 8, path 49, frame 18 image, to illustrate the High-Change Layer (HCL) mask and processing units. The pre-burn Landsat-7 image (a) was paired with the post-burn image (b) to generate differenced Normalized Burn Ratio (dNBR) (c), and the HCL mask (d). The HCL mask was intersected with hotspots (e) and converted to polygon format, to be used as processing units (f) to map the burned area of individual fire events of a selected year (2013).

All subsequent image processing is undertaken within each HCL processing unit, using statistics derived from pixels within the unit. Thus, the HCL mask serves as a sampling region to compute a mean dNBR value, which is then subtracted from the dNBR value of each pixel within the HCL processing unit to create a dNBR differential metric (Key 2005). The burn threshold value is two standard deviations above the mean dNBR differential metric from the sampling region of the HCL mask, which sets the lowest dNBR value for unburned pixels within the HCL processing unit (Fig. 4). Pixel values with a dNBR differential metric less than the burn threshold are classified as burned. This threshold is deemed adaptive, as the resulting mean dNBR and standard deviation are unique for each HCL processing unit from which burned and unburned pixels are identified. Fig. 4 shows an example for HCL processing unit 102 (from Fig. 3f) from which the pixels classified as burned were converted to polygons (Fig. S2, Supplementary material). Burned area polygons are subject to a quality control process that includes visual inspection and post-editing to correct for errors, if any, in the automated delineation. This quality control is achieved using a digitally enhanced Landsat post-fire image as the backdrop in the GIS (geographic information system), to interpret features where burned area mapping is challenging. For example, fires may occur near or within areas that were previously burned, which results in the need to map one burn over another (Fig. S3, Supplementary material). Hotspots (Fig. S3a) must then be used to adjust what was initially mapped (Fig. S3b) to the area that was actually burned (Fig. S3c). Other examples that may require post-editing include burned areas within wetland vegetation (Fig. S4a, Supplementary material), rocky outcrops (Fig. S4b) and forest cutblocks (Fig. S4c).

Fig. 4.  Example histogram of High-Change Layer (HCL) processing unit 102, to illustrate threshold of burned and unburned pixels at 2 s.d. above the mean of the HCL mask for this burn. dNBR, differenced Normalized Burn Ratio.

The third data source used in NBAC is derived from the HANDS algorithm (Fraser et al. 2000) that Natural Resources Canada uses to generate the first national fire maps at the end of each fire season for all fires at least 250 ha in size. The HANDS system started production in 1998 using SPOT VEGETATION 1000-m resolution imagery, which in 2013 was replaced by higher-spatial-resolution Proba-V imagery resampled to 250 m with a corresponding reduction to 50 ha in the minimum size of burned areas mapped by this method.

The HANDS algorithm uses a normalised NIR-SWIR (near infrared-shortwave infrared) index from pre- and post-burn imagery, fire hotspots and land cover to map burned pixels on the basis of adaptive image-differencing thresholds. The output is a raster-based burn mask that is subsequently converted to polygon format for incorporation into the NBAC spatial data warehouse, a database that contains burned areas mapped from the three sources (Fig. 1). Burned areas will generally be overestimated when mapped from coarse-spatial-resolution satellite data (Eva and Lambin 1998; Fraser et al. 2004). To account for this overestimation, area-based calibration models based on Theil–Sen regression (Fernandes and Leblanc 2005) are used to generate statistically unbiased national burn area estimates from HANDS fire polygons for each year of the NBAC. Spatial and temporal filters are also applied to reduce commission errors by removing fires that were less than 250 ha in size, as well as fires with a start date before 16 April or after 30 September (as determined from MODIS hotspots) and fires with burned area consisting of more than 50% agricultural land cover.

In the NBAC spatial data warehouse, each fire event is assigned a unique National Fire Identification (NFIREID) value (Fig. 1), where a fire event comprises one or more polygons over a localised spatial region for a specific year. The unique NFIREID is then assigned to agency, MAFiMS and HANDS fire polygons associated with that fire event. Through this process, an annual fire event layer is derived from a 1250-m buffered radius of the HANDS fire polygons, augmented by buffered agency and MAFiMS polygons for fire events lacking HANDS coverage, as these various sources are incorporated or mapped. Thus, agency, MAFiMS and HANDS polygons intersecting with the same fire event polygon in a given year are assigned the same NFIREID value, which results in a multisource mapping product for a single fire event. Fire records in the attribute table may also contain duplicate NFIREID values where partial burn values are mapped, and where agency fire polygons are split by jurisdictional boundaries in which each agency mapped its region with a different data source or method. In these examples, the different partial burn, agency, source and method values are retained for the same NFIREID.

The NBAC system applies user-defined decision rules considered to be the ‘best source’ to assemble the national product for each fire event (Fig. 1). Therefore, for each NFIREID, the NBAC product will retain only a single mapping source from among those available, according to the selection priority assigned to each source within the rule set. In general, the NBAC selects MAFiMS polygons when available, then agency polygons, followed by HANDS polygons (if no other data source is available). This order of precedence depends on how well the delineation follows the burn pattern observed in the post-fire image and coincident MODIS hotspots within the perimeter, as well as on the data source and method. The rules are, however, flexible. For example, if a particular agency polygon was derived from higher-spatial-resolution data (from the RapidEye and QuickBird sensors) or from aerial photography, then a subrule specifies that these polygons should be preferentially selected. Each NBAC polygon is assigned attributes that either describe the fire itself or that provide metadata about how the event was mapped (Table 2). Burned area, fire cause and fire dates describe features of the fire event, whereas the data source, date of image capture and mapping method are metadata that provide insights into the quality of the mapped fire polygon. Fire start and end dates are derived from both MODIS hotspots and agency data, when available. Fire causes are transferred from agency polygons to MAFiMS fire polygons in the NBAC. If desired, the NBAC can also be used to derive single-source national products with a consistent set of attributes and metadata.

Table 2.  Subset of descriptive fire attributes associated with each fire polygon in the National Burned Area Composite

The accuracy of Landsat-derived burned areas generated from MAFiMS was previously unknown. Undertaking an accuracy assessment requires independent reference data that can be considered more accurate than the product being assessed (Olofsson et al. 2014; Sparks et al. 2015). Agency data considered to be suitable reference data were burned polygons derived from interpretations of aerial photographs from soft-copy photogrammetry and high-spatial-resolution satellite imagery from the RapidEye sensor. These two sources alone, however, would result in only a limited sample of independent polygons. Visual interpretation of Landsat images has also been used as reference data when performed by a well-trained expert (Sparks et al. 2015; Vanderhoof et al. 2017).

A total of 41 burned area pairs from MAFiMS and agency fire polygons were compiled for the assessment. This sample size was low, because MAFiMS polygons are not usually produced when high-quality agency polygons are deemed to accurately and precisely represent the fire event. For these selected fires, the agency polygons were from the provinces of Alberta in western Canada (Boreal Plains Ecozone), and Quebec (Boreal Shield Ecozone) in eastern Canada. The burn pairs represented different vegetation types to the extent these two regions could considering they occurred in different ecological regions

Remote sensing-based burned area products may come with measures of accuracy (Domenikiotis et al. 2002; Henry 2008; Polychronaki and Gitas 2012; Guindon et al. 2014). Further insights into uncertainty could be gained by also estimating measures of precision (confidence intervals) for the derived area burned statistics, but these are rarely calculated. For a primary use of NBAC data, Canada’s NFCMARS, confidence intervals are critical because good practice guidance from the Intergovernmental Panel on Climate Change requires GHG inventories to identify, quantify and reduce uncertainties in reported estimates (Penman et al. 2003). Analysis for the NFCMARS (Metsaranta et al. 2017), as reported in Canada’s GHG inventory, assumes an uncertainty of ±10%, for NBAC-derived area burned statistics, but this value is based on judgement, not analysis. It is therefore important to determine if actual uncertainty is larger or smaller than currently assumed.

Different data sources and methods may give different estimates of area burned for fire events (Parisien et al. 2006), potentially forming the basis for using empirical Monte Carlo simulation approaches for confidence interval estimation. Various data sources and methods in the NBAC were placed into six confidence classes, ranked from lowest to highest uncertainty as follows: Class 1, air photo, field GPS, and satellite RapidEye and QuickBird sensors; Class 2, the MAFiMS; Class 3, agency satellite with 30-m resolution (Yukon, Saskatchewan, Manitoba, Quebec) or 20-m resolution (Alberta (Sentinel-2)); Class 4, airborne GPS (British Columbia, Ontario) and SPOT-4 (Saskatchewan); Class 5, HANDS-calibrated SPOT Vegetation and Proba-V; and Class 6, undefined, buffered hotspots and GPS, sketch maps or MODIS polygons. In the NBAC, 32% of the 16 139 fires were mapped by more than one confidence class method. Few fires (n = 41) using Class 1 were also mapped by MAFiMS (Class 2); therefore, Classes 1 and 2 were merged and deemed the reference class where confidence in fire perimeters was highest. Classes 3 to 6 were compared with this merged reference class. The number of fires mapped by both a given confidence class method and the reference method ranged from 976 for Class 3 to 2593 for Class 5; the sample sizes were intermediate for Class 4 (n = 795) and Class 6 (n = 771).

Confidence intervals for area burned statistics were derived using empirical Monte Carlo simulation methods. The sample of fires mapped by a given confidence class method and the reference method were used to estimate an error distribution for that method. A simple proportion error (Ej) for the estimated area (A) for each fire j was calculated against the reference area A_R for fire j, as follows:

The distribution of error values for all fires mapped by each confidence class (C) was used as an estimate of the overall error distribution for that confidence class (E_C). As expected, errors varied by confidence class overall, with larger errors for the lowest confidence class and smaller errors for the highest confidence class. However, errors were also specific to the size class, tending to be larger for smaller fires (Fig. 5). For the estimation of confidence intervals, both size and method class were considered in the calculation. This calculation involved randomly generating estimates of the error (E) for the estimated area burned (A) by each fire j in the NBAC; this error value was then used to generate an alternative estimate of the area burned that accounted for the observed error for the combination of size and confidence class of fire j, with this process being repeated for 1000 iterations for each fire event. The errors were modelled using the empirical cumulative distribution function of each combination of confidence class and size class. The results for each iteration were summed nationally by year, and approximate 95% confidence intervals were derived from the 2.5th and 97.5th percentiles of the distribution of values obtained from the 1000 simulations. The resulting estimates were compared with the default value used in uncertainty assessments by the NFCMARS (±10%) for GHG inventory reporting (Metsaranta et al. 2017) to determine if they were more or less than what is currently assumed.

Fig. 5.  Distribution of error values for each confidence class fire mapping method used in the National Burned Area Composite (NBAC). Confidence classes are ranked from lowest (Class 1) to highest (Class 6) uncertainty. Class 1, which is omitted from the figure because it has no reference data, includes fires mapped by air photo, field global positioning system (GPS), and satellite RapidEye and QuickBird sensors. Class 2 (panel (a)) includes fires mapped by the Multi-Acquisition Fire Mapping System. Class 3 (panel (b)) includes fires mapped by agency satellite with 30-m resolution (Yukon, Saskatchewan, Manitoba, Quebec) or 20-m resolution (Alberta (Sentinel-2)). Class 4 (panel (c)) includes fires mapped by airborne GPS (British Columbia, Ontario) or SPOT-4 (Saskatchewan). Class 5 (panel (d)) includes fires mapped by Hotspot and Normalized Difference Vegetation Index Differencing Synergy (HANDS) calibrated SPOT Vegetation and European Space Agency Proba-V. Class 6 (panel (e)) contains all other mapping methods, including undefined, buffered hotspots and GPS, sketch maps or Moderate Resolution Imaging Spectrophotometer (MODIS) polygons. For Classes 3 through 6, confidence Classes 1 and 2 are used as the reference data. For Class 2, only Class 1 is used as the reference data. Squares represent the median error value, vertical lines the 5th and 95th percentiles (90% of values), and circles the 2.5th and 97.5th percentiles (95% of values). Within each confidence class, the errors are size-specific, and the values for four size classes, representing fires with areas burned in the 0 to <25%, ≥25 to <50%, ≥50 to <75% and ≥75 to 100% sample percentiles, are plotted individually, from left to right in each panel, using the labels 0–25, 25–50, 50–75, and 75–100 respectively. For clarity, the y-axis is truncated at a maximum of 5, representing an estimated area that is five times the reference area. For confidence Classes 5 and 6, the 97.5th percentile for the smallest size class would be 18 and 39 respectively.

For the period 2004 to 2016, annual area burned statistics in the NBAC vary considerably from year to year (Table 3), consistent with previous reports of high interannual variability in burned areas across Canada (e.g. Stocks et al. 2003; Hanes et al. 2019). Across these 13 years, the total area burned, as reported by the NBAC, ranged from 0.77 Mha in 2009 to 4.2 Mha in 2013, with an average annual area burned of 2.26 Mha. On average, the largest proportion of burned area was mapped by the MAFiMS (65.6%); most of the remaining area was mapped from agency data (32.9%), with very little data mapped from HANDS data (1.5%). In some years, however, half or more of the area burned was accounted for by agency data. The average annual proportion of burned area from all MAFiMS and agency fire polygons in the NBAC that used data with spatial resolution of 30 m or finer was 86.3% (Table 3). Across all burned areas and years, unburned forest, water or both was removed from ~90% of the burned area (representing 17% of the number of polygons) in the NBAC (Table 3). Perimeter-based maps of wildfires that include unburned forest and water overestimate area burned (Kolden and Weisberg 2007) and thus, their removal in NBAC results in more accurate estimates of actual area burned. Visual indications of where fires are occurring in Canada can be obtained from the spatial distribution of NBAC polygons across the boreal and non-boreal zones and in relation to the portion of Canada’s forest considered to be managed for GHG reporting purposes (Stinson et al. 2011; Metsaranta et al. 2017; Kurz et al. 2018). Approximately 96% of the burned area in Canada from 2004 to 2016 occurred in the boreal zone, with more than 50% occurring in the non-managed forest and 45% in the managed forest of this zone (Fig. S5, Supplementary material).

Table 3.  The proportion of burned area in the National Burned Area Composite, by source, generation from high spatial resolution data (≤30-m), and burned area for which unburned forest and/or water have been removed, 2004 to 2016

Lightning was the leading cause of fires in the NBAC in terms of both number of fire events (58.4%) and proportion of burned area (86%) (Fig. S6, Supplementary material), similarly to results reported by Hanes et al. (2019). Human-caused fires accounted for 31.2% of the fire events, but only 2.6% of the total burned area. Fires caused by lightning occurred across all of the boreal forest, but were most prominent in western Canada including the Boreal Shield, Boreal Plains, Taiga Plains and Taiga Shield Ecozones (Fig. S6; Hanes et al. 2019). Human-caused fires occurred in the more southern ecozones and tended to be more scattered in distribution. The undefined category of fire cause, derived from fire polygons in which agencies did not label a cause to the fire, represented ~10% of fire counts and burned area in the NBAC.

The burns sampled ranged in size from 254 to 124 687 ha, with the MAFiMS polygons ranging from 255 to 125 976 ha, and they were highly statistically correlated with the agency polygons for the same fires (r = 0.998, P < 0.001). The absolute value of the difference in burned area between agency and MAFiMS data was used to compute an average accuracy of 96% (s.d. ±3%). Another way to express the relation between agency and MAFiMS burned area is to generate a continuous linear model (agency burned area = 561.7 + 0.991 × MAFiMS burned area); this resulted in a high coefficient of determination (adjusted r² = 0.996), with a computed prediction error sums of squares that resulted in a root-mean-square error of 2 284 ha. A caveat for such a model is that the points in the regression represent fire events with different burned areas, but the relative residuals for events larger than the median event (−0.01) were found to be smaller than the relative residuals for events smaller than the median event (0.22). Therefore, on average, the agreement between the reference (agency) and MAFiMS products increases with the size of the fire event. Given that the analysis included considerably more small fire events than large ones, this correlation estimate is conservative, and we can conclude that the accuracy of MAFiMS is very close to that of the best manual delineations derived by the agencies.

The proportion of fires in each year as mapped by the six confidence classes is shown in Fig. 6a with the resulting confidence intervals for area burned shown in terms of area (ha) in Fig. 6b and as a proportion of annual area burned in Fig. 6c. The confidence intervals shown in Fig. 6b and 6c represent the full span from the 2.5th to the 97.5th percentile of the 1000 Monte Carlo simulations. The total width of the confidence intervals spanned on average 186 565 ha or 8.7% of the annual area burned, varying interannually from a minimum of 3.5% (81 705 ha) to a maximum of 16.6% (494 256 ha). The annual variation in calculated confidence intervals was due to differences in the proportion of fires mapped in each year by different confidence class methods. The width of the confidence interval is related to the proportion of fires mapped using the confidence Class 1 and 2 methods, with the width being narrower when more area is mapped using the high confidence class methods (r² = 0.35, P < 0.05). Fig. 6b and 6c also show confidence interval estimates for each year based on the default assumption currently used for uncertainty analysis for GHG inventory reporting, ±10% (Metsaranta et al. 2017). The calculated confidence intervals were on average ±4.3% of the annual area burned in a given year, which was approximately half the default assumption. There was considerable variation in the numerical values for the estimated 95% confidence intervals for annual areas burned (Table 3).

Fig. 6.  Proportion of area in the National Burned Area Composite (NBAC), mapped by each confidence class method in each year (a) and confidence intervals for annual areas burned, as calculated in the present study (b and c). The values plotted in panels (b) and (c) represent the total width spanned by the 2.5th and 97.5th percentiles of 1000 Monte Carlo simulations accounting for confidence and size class errors for each individual fire in the NBAC. Panel (b) shows confidence intervals in area terms, panel (c) as a proportion of annual area burned in a given year. The legend next to panel (a) applies only to panel (a). The legend next to panel (c), which applies to both panels (b) and (c), indicates the confidence intervals calculated in the present study as dark grey bars and labelled ‘Calculated’, and the confidence intervals currently assumed in uncertainty analyses for Canada’s National Forest Carbon Monitoring, Accounting, and Reporting System (Metsaranta et al. 2017) are shown as light grey bars and labelled ‘Default’.

The present paper has described the NBAC, which generates an annual fire disturbance product for input to the NFCMARS (Stinson et al. 2011; Metsaranta et al. 2017). Three questions relevant to the NBAC merit discussion: (i) what are the advantages of the NBAC? (ii) What is the uncertainty associated with estimates of area burned? (iii) How can the NBAC be updated and its time series extended?

Fire management agencies report the area burned in polygon format (Stocks et al. 2003; Hanes et al. 2019). Therefore, during development of the NBAC, it was decided to represent the spatial unit of burned area as a polygon to facilitate integration with agency data. With all of the NBAC data sources being sequentially built as polygon datasets, the NBAC product becomes capable of storing descriptive information about fire events (e.g. fire cause, fire start and end dates, data source, mapping method) (Table 2). Providing the burned area in polygon format, with a set of attributes, enables spatial database queries to extract fire information by specific attributes that could support a multitude of ecological studies. In addition, the NBAC system was made adaptable to incorporate new or updated information from various data sources that could be added to the NBAC attribute table. The application of polygon features in NBAC also facilitates manual editing of the MAFiMS perimeter boundary in a GIS if improvements in delineation around the burn scar are required (Fig. S4).

A key concept in the development of the NBAC was the selection of the ‘best data source’ or any user-defined source prioritisation, from multisource to single source, for representation of a burned area product. Using only remote sensing data with 30-m spatial resolution from satellites such as the Landsat Enhanced Thematic Tracker (ETM+) or Operational Line Imager (OLI) results in consistency of the data source, but satellite imagery is not necessarily the best data source. We found that generation of agency burn polygons by soft-copy photogrammetric workstations using digital aerial photographs and high-spatial-resolution satellite data (such as data from RapidEye and QuickBird sensors) yielded results that were as precise or more refined than could be derived from 30-m satellite data. As a result, there was often no justifiable basis to exclude agency burned area data derived from high-spatial-resolution data sources during creation of the NBAC (Table 3).

Mapping based on remote sensing data can be challenging when applied over a diverse range of fire-affected environments (Lentile et al. 2006). The state of burned vegetation and its subsequent spectral response depend on the pre-fire fuel characteristics and burning (or fire weather) conditions that occurred during the fire (de Groot et al. 2007). The spectral and spatial heterogeneity of burned patches can complicate automatic detection algorithms applied to satellite images (Chu and Guo 2014). The image mapping approach developed for the NBAC largely addresses the variability of burned area within a fire by using the automated iterative, adaptive threshold approach in MAFiMS. Certain cover types, however, can confound the spectral signature of burned pixels from Landsat imagery, which may necessitate post-editing portions of the polygon boundary (Fig. S4). Selecting the acquisition date for a post-burn image can pose challenges depending on when fire events start and end due to the spectral responses from post-burn vegetation recovery. An example from Landsat OLI imagery is a fire that burned from 5 to 9 May 2016 is clearly visible in an image acquired on 17 May 2016 (Fig. S7a, Supplementary material). By mid-summer (20 July 2016), the burn scar had a diminished, dark reddish tone (Fig. S7b), and by 22 September 2016, the greenish tones associated with vegetation recovery appeared similar to the surrounding forest vegetation (Fig. S7c). MAFiMS employs fire hotspots for selection of post-fire imagery to allow determination of the appropriate fire start and end dates. From meeting these kinds of challenges, could the NBAC be considered the data source for national burned area reporting?

The CNFDB compiles annual fire data that includes both polygon and point data from the agencies, which are completed several months after the fire season. Since the first production of the NBAC (for the 2004 fire season), its estimates of area burned have been mostly lower than those of the CNFDB, largely because the MAFiMS was used to replace large, broadly delineated agency polygons that included unburned forest islands and water. Only for a few years, such as 2007, was the burned area estimated by NBAC larger than that estimated by the CNFDB, because the NBAC captured fire events that were not mapped by the agencies. Provincial fire management agencies also do not map all fires with polygons in every year that will contribute to the differences observed (Table S1). The burned areas reported in CNFDB statistics and the Canadian Interagency Forest Fire Centre (CIFFC) Canada Report (https://ciffc.ca/publications/canada-reports accessed 30 July 2020) are not derived from the final post-edited agency polygons or from MAFiMS fire products, and therefore have tended to be larger than those reported in the NBAC (Fig. 7).

Fig. 7.  Annual burned area statistics from 2004 to 2016 derived from various sources: National Burned Area Composite (NBAC; present paper), National Forestry Database (NFD; Natural Resources Canada 2019c), Canadian National Fire Database (CNFDB; Natural Resources Canada 2019d) and Canadian Interagency Forest Fire Centre (CIFFC; https://www.ciffc.ca). The error bars for the NBAC are 95% confidence intervals generated by empirical Monte Carlo simulation methods, as described in the present paper. Confidence interval estimates for the other sources of national area burned statistics have never been calculated, so they are unknown.

Having several sources for compiling national-level fire data summary reports and statistics can be problematic, particularly in defining the burned area in Canada. In recognising the problem, an initiative is under way through CIFFC partner agencies to establish national data standards and to develop a national information management strategy (CIFFC 2015). It has been recognised through this initiative that the NBAC may be well suited to be integrated into this pan-Canadian data integration and governance strategy.

Systems for monitoring wildfire, including the NBAC, are associated with uncertainties attributable to various systematic and random errors that influence both the accuracy and precision of the derived area burned statistics. Error sources have been discussed in previous national-scale fire mapping studies in Canada (e.g. Stocks et al. 2003; White et al. 2017; Hanes et al. 2019), but they have never before been formally quantified. Thus, the calculation of numerical confidence bounds for national burned area statistics derived from the NBAC was a unique contribution of the current study. Alternative sources of national-scale area burned statistics (the National Forestry Database (Natural Resources Canada 2019c), the CNFDB (Hanes et al. 2019; Natural Resources Canada 2019d) and the CIFFC (https://www.ciffc.ca/ accessed 20 July 2020)) do not provide these estimates. As a result, although alternative sources of statistics can have different values in various years, inferring the cause of differences and their significance remains speculative because confidence interval estimates are lacking for sources other than the NBAC.

The approach used for calculating confidence intervals for NBAC statistics was highly empirical, because it used Monte Carlo simulation and the observed difference in estimated areas between different methods of mapping individual fire events of different size and confidence class as the basis for the propagation of errors. The magnitude of estimated uncertainties depends on the data and assumptions used. Alternative approaches might yield different results, and therefore, true uncertainty may differ from what is reported here. For example, the analytical process considered estimated areas only, not differences in the alignment of burned area perimeters. This is one factor that could cause estimated uncertainty to differ from true uncertainty. Although no true data providing the exact area burned exist, we used the best data available (Congalton and Green 2008). In the present study, 41 fires were available for the best reference class, against which estimates from the MAFiMS could be compared. Thus, the MAFiMS itself was used as the reference class for judging fires mapped using lower-confidence methods. Increasing the size and representativeness of this reference set would improve the estimates of accuracy and uncertainty (Hawbaker et al. 2017; Vanderhoof et al. 2017). This would not alter the conclusion that NBAC-derived area burned statistics are fairly precise. For nearly all years, calculated confidence intervals were narrower (±4.3%) than the ±10% uncertainty assumed for GHG inventory reporting (Metsaranta et al. 2017). In other words, relative to uncertainty about other ecosystem processes such as forest productivity and soil carbon dynamics, uncertainty about area burned by wildfire contributes relatively little uncertainty to current estimates of the carbon balance of Canada’s forests (Metsaranta et al. 2017). Incorporating improved methods, particularly in measuring variation in fire severity, will further reduce this uncertainty, but for the purposes of GHG reporting, both area burned and its uncertainty can be considered well characterised. Metsaranta et al. (2017) provides further details regarding the contribution of various factors to uncertainty in forest GHG inventory reporting for Canada.

The NBAC system can be adapted to new image sensor opportunities and fire data sources. In 2013, the MAFiMS was updated to allow preprocessing of Landsat OLI imagery for pairing with Landsat Thematic Mapper (TM) and ETM+ gap-filled imagery for burned area mapping. MAFiMS fire polygons derived from Sentinel-2 imagery were recently implemented for the 2019 fire season and work is currently under way to process VIIRS (Visible Infrared Imaging Radiometer Suite) fire hotspots for MAFiMS application. The transition from SPOT VGT (Vegetation) to Proba-V imagery for HANDS has also improved the data quality of the few fire events mapped by HANDS that remain in NBAC, owing to the lack of delineations from higher-spatial-resolution data. Some agencies also continue to improve their data quality by updating historical fire polygons derived from aerial surveys to new fire polygons using Landsat imagery. With each release of updated agency fire polygons, the data are replaced in the NBAC system to improve national statistics. The Canadian Forest Service recently partnered with the Government of Northwest Territories to update over 30 years of fire polygons largely derived from aerial GPS tracking. These additional MAFiMS fire polygons of Northwest Territories fires have been incorporated into the NBAC, along with other newly mapped MAFiMS fire polygons of large fires mapped across the country and modified data from the National Terrestrial Ecosystem Monitoring System Composite2Change project (White et al. 2017), which extends the NBAC time series to 1986. This extended time series results in a data record of more than 30 years of the best available fire polygon data that can help support science, policy and national reporting (Sankey 2018). Current versions of the NBAC are freely available on the CWFIS website (http://cwfis.cfs.nrcan.gc.ca/datamart accessed 20 July 2020).

There has been considerable interest in mapping burned areas caused by forest fires at regional and national levels. The many data sources and methods employed have resulted in considerable variation in terms of delineating particular burned areas and the uncertainty associated with each estimate. In Canada, fire management agencies may employ more than one data source and method within their respective jurisdictions, and there are differences in the degree to which unburned forest islands and water features are removed from within the burn perimeter. To address these inconsistencies, the NBAC was developed to assign the best data source for each given burn event. This system has been functioning on a continuing annual basis since its inception in 2004, and for the purposes of this paper, data were extracted across a 13-year time series (from 2004 to 2016).

Mapping burned areas as polygons has enabled integration with agency data, from which a set of descriptor attributes is also available for each burn event contained within the national product. A key concept associated with compiling polygons of the same fire event was the opportunity, during creation of the NBAC, to select from among multiple data sources in the mapping of fires in the year of occurrence. This concept sets the NBAC apart as a unique national burned area dataset, distinct from pixel-based disturbance products.

Across the 13-year time series, ~90% of the burned area in the NBAC is represented by polygons where unburned forest islands, water features, or both, are removed. A process for assigning levels of confidence to the data source and mapping method was used to derive a measure of uncertainty for each year of the NBAC. Results from this process demonstrated that calculated uncertainty was consistently smaller (±4.3%) than the previous fixed assumption of ±10% for the NBAC.

The authors declare no conflicts of interest.