Four National Maps of Broad Forest Type Provide Inconsistent Answers to the Question of What Burns in Canada

Wildfires are burning increasingly extensive areas of forest in Canada, reducing their capacity as carbon sinks. Here we compare the answers that four independent land cover datasets, produced from different satellite images (SPOT, Landsat, and MODIS), provide for the question of what burned in Canada in recent years. We harmonized the different datasets into a common, simpler legend consisting of three classes of forest (needle-leaf, broadleaf, and mixed) plus non-forest, and resampled them to a common pixel size (250 m). Then we used annual maps of burned area to count, for each map and year from 2011 to 2014, the number of burned pixels of each class, and we tabulated them by terrestrial ecozone and Canada-wide. While all four maps agree that needle-leaf forest is the most frequently burned class in Canada, there is great disparity in the results from each map regarding the proportion of burned area that each class represents. Proportions reported by one map can be more than double those reported by another map, and more than four times at the ecozone level. We discuss the various factors that can explain the observed discrepancies and conclude that none of the maps provides a sufficiently accurate answer for applications such as carbon accounting. There is a need for better information in areas lacking forest inventory, especially in the vast unmanaged forest of Canada.


Introduction
Article 5 of the Paris Agreement requires signatory countries to conserve and enhance sinks and reservoirs of greenhouse gases (GHG), including forests [1].Wildfires reduce the capacity of forests to sequester carbon, to a point where the combined effects of insects and fire can transform a forest region from a carbon sink to a source [2].This is especially relevant in the boreal region of Canada, where fires burn an average of over two million hectares (Mha) a year [3], and where climate change is expected to increase the frequency, severity, and extent of both pest outbreaks and wildfires [4,5].Reliable estimation of GHG emissions and other losses caused by fire requires spatially explicit, accurate data on forest attributes from which to derive fuel type and load [6,7].Such data are only available, however, for portions of the 230 Mha of Canada's forests that are managed, whereas for the 120 Mha that remain unmanaged, estimates mostly rely on land cover information derived from satellite imagery [8].Here we use four such data sources, plus a Canada-wide annual map of burned area, to derive estimates of burn area by four broad categories, namely non-forest (e.g., shrubland, grassland), needle-leaf forest, broadleaf forest, and mixed forest.The objective of this study was to compare the estimates from each data source to determine whether the results regarding what burned in Canada in a given year were consistent across the sources.Several papers in the literature have compared large area land cover maps (e.g., [9][10][11]), but to the best of our knowledge, this is the first where the comparison deals with the estimation of the distribution of burned area across different forest types.This latter aspect has indeed been studied before, but using a single data source (e.g., [12]).

Materials and Methods
The 6.3 million km 2 (626 Mha) study region comprises the nine terrestrial ecozones of Canada [13] that are mostly covered by forest; thus it excludes arctic and prairie ecozones (Figure 1).It also excludes some half million km 2 that lacked data in some of the sources, mainly in the Atlantic Maritime ecozone in Nova Scotia and Prince Edward Island plus the Mixedwood Plains ecozone in southern Ontario; these are areas that in any case are rarely affected by fire.
in Canada in a given year were consistent across the sources.Several papers in the literature have compared large area land cover maps (e.g., [9][10][11]), but to the best of our knowledge, this is the first where the comparison deals with the estimation of the distribution of burned area across different forest types.This latter aspect has indeed been studied before, but using a single data source (e.g., [12]).

Materials and Methods
The 6.3 million km 2 (626 Mha) study region comprises the nine terrestrial ecozones of Canada [13] that are mostly covered by forest; thus it excludes arctic and prairie ecozones (Figure 1).It also excludes some half million km 2 that lacked data in some of the sources, mainly in the Atlantic Maritime ecozone in Nova Scotia and Prince Edward Island plus the Mixedwood Plains ecozone in southern Ontario; these are areas that in any case are rarely affected by fire.

Burned Area Data
We used the four polygon vector layers of the National Burned Area Composite (NBAC) that correspond to each of the years 2011-2014 respectively, totaling 11.6 Mha (yellow blobs in Figure 2).NBAC provides delineations of areas burned by fire events >10 ha across Canada in each year since 2004 and uses the best available information for each fire event, which typically is based on a semiautomated procedure that uses pre-and post-fire Landsat imagery [8,14].Each vector layer was converted to a 30 m pixel raster mask in the NAD1983 Albers equal-area conic projection.While the period chosen may not be representative of what would burn in Canada over a longer period, we wanted to ensure that all land cover data sources reflected the pre-fire cover conditions (two of the four sources represent conditions as recent as 2010).

Burned Area Data
We used the four polygon vector layers of the National Burned Area Composite (NBAC) that correspond to each of the years 2011-2014 respectively, totaling 11.6 Mha (yellow blobs in Figure 2).NBAC provides delineations of areas burned by fire events >10 ha across Canada in each year since 2004 and uses the best available information for each fire event, which typically is based on a semi-automated procedure that uses pre-and post-fire Landsat imagery [8,14].Each vector layer was converted to a 30 m pixel raster mask in the NAD1983 Albers equal-area conic projection.While the period chosen may not be representative of what would burn in Canada over a longer period, we wanted to ensure that all land cover data sources reflected the pre-fire cover conditions (two of the four sources represent conditions as recent as 2010).

Land Cover Data
We obtained, in raster format, the following four land cover datasets from their authors or from publicly accessible repositories.We briefly report on source input imagery, acquisition period, classification method, number of classes, and thematic accuracy for each dataset.

EOSD Land Cover
The Earth Observation for Sustainable Development (EOSD) 25 m pixel raster is a circa 2000 national land cover map of the forested area of Canada [15].The project was initiated as a partnership between the Canadian Forest Service and the Canadian Space Agency, with support from provincial and territorial governments.The product was derived from digital classification of Landsat-7 ETM+ imagery acquired mostly in the summer of 2000, but in cloudy areas the image could be from up to two years earlier or later.Spectral clusters were generated through the k-means algorithm and then merged manually into similar cluster groups, where each group was assigned to one of the 21 land cover classes in the EOSD legend.Although the EOSD product has not been rigorously validated, a regional accuracy assessment was undertaken in Vancouver Island [16].Airborne video data were acquired over 15,000 ha along four flight lines parallel to the main axis of the island.The center of 500 video frames randomly selected from a systematic sample of 2500 frames spaced 15 seconds apart was photo-interpreted in a 90 m square around the center, and the interpreters' calls were compared with the majority EOSD label of pixels inside that square [16].The reported overall accuracy when tree density in the forest classes was not considered (i.e., 15-class legend) was 73%.

Medium-Resolution Land Cover Mapping of Canada from SPOT 4/5 Data
This 20 m land cover map of Canada's forested regions was generated from SPOT 4/5 imagery by the Canada Centre for Mapping and Earth Observation [17].Over 5000 SPOT 4-5 summer images (60 by 60 km) acquired from 2005 to 2010 were collected and preprocessed to balance image radiometry to create 20 mosaicked mapping zones that were classified separately using a common 16-class legend.Classification consisted of a fuzzy k-means clustering algorithm followed by manual merging and labeling of clusters.The accuracy assessment was based on 1566 reference points collected from a purposive field survey and from random sampling and visual interpretation of submetric Google Earth imagery.The reported overall accuracy was 71%.

Canada's MODIS Land Cover Time Series (LCTS)
This series of annual (2000-2011) land cover maps of 250 m pixel size was developed by the Canada Centre for Mapping and Earth Observation using satellite data from the MODIS sensor [18].

Land Cover Data
We obtained, in raster format, the following four land cover datasets from their authors or from publicly accessible repositories.We briefly report on source input imagery, acquisition period, classification method, number of classes, and thematic accuracy for each dataset.

EOSD Land Cover
The Earth Observation for Sustainable Development (EOSD) 25 m pixel raster is a circa 2000 national land cover map of the forested area of Canada [15].The project was initiated as a partnership between the Canadian Forest Service and the Canadian Space Agency, with support from provincial and territorial governments.The product was derived from digital classification of Landsat-7 ETM+ imagery acquired mostly in the summer of 2000, but in cloudy areas the image could be from up to two years earlier or later.Spectral clusters were generated through the k-means algorithm and then merged manually into similar cluster groups, where each group was assigned to one of the 21 land cover classes in the EOSD legend.Although the EOSD product has not been rigorously validated, a regional accuracy assessment was undertaken in Vancouver Island [16].Airborne video data were acquired over 15,000 ha along four flight lines parallel to the main axis of the island.The center of 500 video frames randomly selected from a systematic sample of 2500 frames spaced 15 seconds apart was photo-interpreted in a 90 m square around the center, and the interpreters' calls were compared with the majority EOSD label of pixels inside that square [16].The reported overall accuracy when tree density in the forest classes was not considered (i.e., 15-class legend) was 73%.

Medium-Resolution Land Cover Mapping of Canada from SPOT 4/5 Data
This 20 m land cover map of Canada's forested regions was generated from SPOT 4/5 imagery by the Canada Centre for Mapping and Earth Observation [17].Over 5000 SPOT 4-5 summer images (60 by 60 km) acquired from 2005 to 2010 were collected and preprocessed to balance image radiometry to create 20 mosaicked mapping zones that were classified separately using a common 16-class legend.Classification consisted of a fuzzy k-means clustering algorithm followed by manual merging and labeling of clusters.The accuracy assessment was based on 1566 reference points collected from a purposive field survey and from random sampling and visual interpretation of submetric Google Earth imagery.The reported overall accuracy was 71%.

Canada's MODIS Land Cover Time Series (LCTS)
This series of annual (2000-2011) land cover maps of 250 m pixel size was developed by the Canada Centre for Mapping and Earth Observation using satellite data from the MODIS sensor [18].Note that for this study, we used the 2010 map.The MODIS LCTS was generated with a backdating/updating approach using annual change detection from a baseline 2005 MODIS land cover map of Canada that was in turn derived from a decision tree classifier applied to monthly composites of selected MODIS bands acquired during that year [19].This baseline map contained 39 classes, and it was reduced to 25 and 19 class versions for development of the MODIS LCTS.Areas of detected change were differentiated as negative or positive, referring to vegetation loss, such as from harvesting and fire, or to vegetation establishment and succession, respectively.The new class assigned to each area of change was that with the maximum support from combined evidence.The latter came from three sources: class memberships generated from a decision tree classifier, class temporal transition likelihoods, and regional and neighborhood class proportions.To assess the accuracy of land cover time series, three approaches were implemented: visual assessment, comparison with a temporally pooled reference sample, and comparison with the MODIS Global LCTS.The comparison with the reference sample resulted in an overall accuracy of 70% at the 19-class thematic resolution.

Attributes of Canada's Forests through k-NN and MODIS
In 2014, the Canadian Forest Service created a series of 250 m rasters of forest attributes for the reference year of 2001.They are based on data from 8202 National Forest Inventory 2 by 2 km photo plots with complete data, which are centered on the nodes of a 20-km grid encompassing Canada's forested area [20].The four two MODIS pixels corresponding to the inside four corners of each photo plot were selected as a reference dataset and forest attributes were extracted, as well as covariates from MODIS spectral data acquired in 2001 and climatic and topographic variables.This reference dataset was randomly split into 70% training and 30% validation subsets.Then the rasters of forest attributes were populated using k-Nearest Neighbour (k-NN) imputation, wherein the six closest training pixels (closest according to the Euclidean distance in the covariate space) are used to impute a value in each empty pixel [20].On the basis of the validation pixels, the percent treed, percent needle-leaf, and percent broadleaf products generated resulted in R 2 values between 0.59 and 0.62, with RMSE values between 25.5 and 18.75 [20].Given that these are continuous attributes, it is not possible to infer the thematic accuracy of the broad land cover classification we derived from these rasters using these error metrics (but see Discussion).

Data Conditioning
First, we harmonized the different legends of each map into a common, simpler legend consisting of only four classes: non-forest (NF), needle-leaf (NL), broadleaf (BL), and mixed forest (MF) (defined in Table 1).This was straightforward except in two instances.The "young forest" class in the SPOT land cover product (about 130,000 km 2 , or 2% of the study region) was assigned to one of the other six forest classes on a per-case basis, depending on the surroundings of the young forest patches.That is, for each individual patch (i.e., set of connected pixels) of young forest, we assigned the forest class that was most frequently found in a one-pixel buffer around the patch.The other instance is the k-NN map, which was based on three k-NN rasters: %treed (% area of the pixel covered by trees), %needleleaf, and %broadleaf.Non-forest class was defined as pixels where %treed <10, while for forested pixels, broadleaf, needle-leaf, and mixed forests were defined respectively as pixels where %needleleaf <25%, >75% and within the interval 25%-75% (same threshold as in EOSD and MODIS LCTS).We note that the EOSD class "wetland treed" was assigned to "needle-leaf," as most treed wetlands in Canada are dominated by black spruce.Second, we resampled the EOSD and SPOT products to 250 m pixels using the majority rule.Third, we intersected each NBAC layer with the layer of ecozones of Canada and converted the result to a 250 m raster using the numeric ecozone ID of intersect polygons as the digital number of the corresponding polygons.Finally, to facilitate comparison, we reprojected all rasters to a common cartographic projection (NAD1983 Albers equal-area conic), making sure that the ouput rasters have all the same dimensions and "snap" to a common grid (i.e., that the NW corner of all rasters is in the exact same location).We also set as "no-data" any pixels for which some of the four sources had no data.To avoid confusion with the original products, we hereafter call these four-class rasters the EOSD, SPOT, LCTS, and k-NN maps, respectively.We acknowledge that these derived maps do not necessarily reflect the quality of the original ones, for different choices in the reclassification or the resampling could have led to different outcomes, plus the new maps contain less thematic information (and in the case of the EOSD and SPOT, less spatial detail).Therefore, our conclusions apply more to our maps than to the original products, one of which (the k-NN) is not even a land cover classification.

Data Analysis
For each map and year 2011-2014, we counted the number of burned pixels within each of the four broad land cover classes and tabulated them nationally and by ecozone.To visually appreciate the differences among answers by year and source to the question of what burned in Canada, we created 16 pie charts (Figure 3), each representing the breakdown of burned area in a given year by cover type according to one of the four maps.To see if the answers differ by ecozone, we aggregated the burned area in those four years and created one bar chart per map where individual bars represent each ecozone and the four segments of each bar represent the percentage of the total burned area in that ecozone that belongs to each of the four land cover classes (Figure 4).To better understand the differences in the answers, we performed a series of additional analyses.First, to have an indication of the contribution of the age difference between maps, we computed the breakdown of area burned in 2011-2014 by cover type according to the first map in the MODIS LCTS (year 2000) in addition to the one we chose (year 2010).Second, to ascertain whether the resampling to 250 m introduced some bias into the proportion of classes in the EOSD and SPOT maps, we computed those proportions before and after the resampling.Third, we assessed whether some inconsistencies regarding non-forest areas could be attributed to different definitions of "forest" in terms of minimum tree cover.For this we computed, using the original k-NN %treed raster, the proportion of pixels in disagreement that would be classified as forest in both the k-NN and the LCTS maps if they had the same definition.Likewise, this time using the %needle-leaf and the %broadleaf rasters in addition to the %treed raster, we computed how many pixels from those that are mixed forest in the SPOT map and either needle-leaf or broadleaf in the k-NN map would become mixed forest if both maps used the same definition.Fourth, to spatially assess disagreements among maps, we also created a raster of agreement among the four maps, where the digital number (DN) of each pixel is the number of maps that assigned the same class to that pixel, ranging from 1 (no agreement, i.e., each map has a different class in that pixel) to 4 (full agreement, i.e., the pixel has the same class in all four maps) (Figure 2).To study how spatial heterogeneity affects agreement, we created a heterogeneity raster per map where the DN of each pixel is the number of pixels in its eight-neighborhood that have a different class; then we averaged the four heterogeneity rasters into a single raster and linearly regressed this combined heterogeneity raster with the raster of agreement.Fifth, to appreciate how the classes in one map get confused with others in a second map, we produced six four by four confusion matrices between each possible pair of maps.We did this four times: (i) using the entire population of pixels in the maps; (ii) using only those pixels inside the set of all burned areas in 2011-2014; (iii) using all pixels inside the managed forest area; and (iv) using all pixels inside the unmanaged forest area.Finally, to ascertain whether the reported accuracy of the original land cover datasets is consistent with the level of agreement we found, we estimated overall accuracy with our four land cover classes for those datasets for which there is a published confusion matrix with the original classes.
Remote Sens. 2016, 8, 539 6 of 13 all four maps) (Figure 2).To study how spatial heterogeneity affects agreement, we created a heterogeneity raster per map where the DN of each pixel is the number of pixels in its eightneighborhood that have a different class; then we averaged the four heterogeneity rasters into a single raster and linearly regressed this combined heterogeneity raster with the raster of agreement.Fifth, to appreciate how the classes in one map get confused with others in a second map, we produced six four by four confusion matrices between each possible pair of maps.We did this four times: (i) using the entire population of pixels in the maps; (ii) using only those pixels inside the set of all burned areas in 2011-2014; (iii) using all pixels inside the managed forest area; and (iv) using all pixels inside the unmanaged forest area.Finally, to ascertain whether the reported accuracy of the original land cover datasets is consistent with the level of agreement we found, we estimated overall accuracy with our four land cover classes for those datasets for which there is a published confusion matrix with the original classes.

Results and Discussion
All four maps agree that needle-leaf forest is by far the forest type most affected by fire in Canada in all years studied and across all ecozones where fire is a common disturbance.Across Canada, the proportion of burned area belonging to this class ranges from 59% (SPOT map, year 2013) to 81% (LCTS map, year 2014) (Figure 3).In all cases, this proportion is greater than what would correspond to a pattern where fire was proportionally distributed among the four classes, meaning that fire has a preference for this class.This is evidenced by the Jacobs selectivity index (JSI, [21]), which ranged between 0.21 and 0.71 depending on source and year (JSI ranges from −1, complete avoidance, to +1, complete preference).While the four maps agree on this issue, any other question raises discrepancies.For example, according to the EOSD map, 36% of burned area in 2013 was non-forest, whereas according to the k-NN map, it was only 13% (Figure 3).Conversely, according to the k-NN map, 12% of burned area in 2014 was mixed forest, whereas according to the EOSD map, it was only 2%.Or take the other two maps: according to the LCTS map, 81% of area burned in 2014 was needleleaf, while SPOT shows 60%; also in 2014, SPOT indicates that 33% of burned area was non-forest, while LCTS shows only 15%.From inspection of Figure 3, several general patterns can be appreciated.For example, k-NN and LCTS have a lower proportion of non-forest in the burned area in all years studied than the other two maps, and the opposite can be said about mixed forest.These general patterns can be explained by the different proportions of each class that the maps exhibit (pie charts in Figure 1).However, there are deviations from the general pattern that defy expectations.For example, in 2011, the difference in mixed forest between the EOSD and LCTS maps is just three percentage points (11% vs. 14%), even though EOSD has just over half the amount of mixed forest that LCTS has (9% vs. 16%).
The discrepancies worsen at the ecozone level even after the area burned in different years is aggregated (Figure 4).For example, according to EOSD, 22% of the area burned in 2011-2014 in the

Results and Discussion
All four maps agree that needle-leaf forest is by far the forest type most affected by fire in Canada in all years studied and across all ecozones where fire is a common disturbance.Across Canada, the proportion of burned area belonging to this class ranges from 59% (SPOT map, year 2013) to 81% (LCTS map, year 2014) (Figure 3).In all cases, this proportion is greater than what would correspond to a pattern where fire was proportionally distributed among the four classes, meaning that fire has a preference for this class.This is evidenced by the Jacobs selectivity index (JSI, [21]), which ranged between 0.21 and 0.71 depending on source and year (JSI ranges from ´1, complete avoidance, to +1, complete preference).While the four maps agree on this issue, any other question raises discrepancies.For example, according to the EOSD map, 36% of burned area in 2013 was non-forest, whereas according to the k-NN map, it was only 13% (Figure 3).Conversely, according to the k-NN map, 12% of burned area in 2014 was mixed forest, whereas according to the EOSD map, it was only 2%.Or take the other two maps: according to the LCTS map, 81% of area burned in 2014 was needle-leaf, while SPOT shows 60%; also in 2014, SPOT indicates that 33% of burned area was non-forest, while LCTS shows only 15%.From inspection of Figure 3, several general patterns can be appreciated.For example, k-NN and LCTS have a lower proportion of non-forest in the burned area in all years studied than the other two maps, and the opposite can be said about mixed forest.These general patterns can be explained by the different proportions of each class that the maps exhibit (pie charts in Figure 1).However, there are deviations from the general pattern that defy expectations.For example, in 2011, the difference in mixed forest between the EOSD and LCTS maps is just three percentage points (11% vs. 14%), even though EOSD has just over half the amount of mixed forest that LCTS has (9% vs. 16%).
The discrepancies worsen at the ecozone level even after the area burned in different years is aggregated (Figure 4).For example, according to EOSD, 22% of the area burned in 2011-2014 in the Boreal Shield (the largest ecozone) was non-forest, whereas the k-NN map says it was just 4%.According to LCTS, 12% of burned area in this same ecozone was mixed forest, whereas the SPOT map says it was just 2% (Figure 4).Among other factors, these regional disparities are a consequence of the long-known fact that the likelihood of misclassification varies markedly from one place to another [22].
Regardless of the two levels of analysis (national and by ecozone), these discrepancies can be explained by a number of factors, namely (1) actual changes on the ground in the period between the image acquisitions of different maps; (2) different definitions of the same class; (3) collocated pixels having different actual footprints on the ground in heterogeneous areas; and (4) classification errors.In the following paragraphs, we speculate on the relative importance of each of these factors.
The different time periods of reference of each map (e.g., from the c. 2000 of EOSD to 2010 for LCTS) could contribute a priori to some of the observed differences, in particular in areas that suffered a stand-replacing disturbance, for which an older map would display a forest class that in a newer map would be non-forest.This would happen, for example, for areas burned in 2003-2009, which would display the pre-fire cover in the EOSD and k-NN maps and the post-fire cover in the LCTS map (note that the SPOT map image acquisition spans 2005-2010, so whether a particular burn area shows pre-or post-fire depends on the image date for that area).However, all areas burned in 2011-2014 will have the pre-fire cover in all maps.It is possible, though, that some of these areas had suffered a previous fire in 2003-2009 and thus could belong to a forest class in the EOSD and k-NN maps and to non-forest in the other maps.This is a remote possibility, however, because re-burn in the same decade is unlikely (according to NBAC, roughly 1% of burned area across Canada in 2004-2014 burned more than once in that period).Therefore, the bulk of true changes (during 2003-2009) from forest to non-forest inside the studied burned areas will mostly be due to deforestation and harvesting.Both factors are practically non-existent in the unmanaged forest and overall, they are responsible for less than 1% of the forest to non-forest conversion per year in the managed forest.Therefore, their contribution to the observed differences can only be negligible.Changes from non-forest to forest could be more common, because it can take as little as 5-10 years for a disturbed area that lost its trees to regain a 10% cover of trees >5 m tall and thus regain its status as forest [23].Therefore, some areas that show as non-forest in the EOSD and k-NN maps because of a fire in the 1990s could indeed appear as forest in the other maps.However, once again because of the rarity of re-burn, this situation should not be frequent inside the 2011-2014 NBAC polygons.Finally, changes in tree species composition leading to change of forest class (e.g., from mixed to conifer) can occur in periods shorter than a decade as a consequence of pests, disease, or dieback if the agent only affects one of the two types of species, but generally this will occur in gaps too small to warrant a change of label in a 250 m pixel.In short, it is unlikely that the older vintage of the EOSD and k-NN maps, combined with the true land cover changes that occurred since they were created, can explain but a slim portion of the differences observed within the burned areas we studied.In that case we could quantify (proportion of each cover type in the area burned in 2011-2014, according to LCTS 2000 vs. LCTS 2010), only the non-forest class changed its proportion by more than one percentage point (13.7% in 2000 vs. 14.9% in 2010).That is, out of the 11.9 Mha that burned in 2011-2014, only 0.14 Mha correspond to areas that suffered a previous change from forest to non-forest in 2000-2010.
Disparities in the definition of the classes by different maps could also be a source of discrepancy.For example, non-forest in EOSD and k-NN is anything with less than 10% tree cover, while in SPOT the threshold is 15% and in LCTS it is 20% (Table 1).The definition of mixed forest in EOSD demands that neither needle-leaf nor broadleaf trees account for 75% or more of total basal area, which is similar to the definition of mixed forest in both k-NN and LCTS LCTS (the three use the same threshold but a different attribute in their original legends: tree cover in LCTS, live above-ground biomass in k-NN, and basal area in EOSD).In contrast, the original SPOT map legend [17] indicates that mixed forest occurs whenever the tree cover of both needle-leaf and broadleaf trees exceeds 2.5% (providing their sum is greater than 15%).These different definitions definitely have an impact.For example, 48% of pixels that are forest in k-NN and non-forest in LCTS fulfill the expectation of having less than 20% tree cover according to the k-NN %treed raster.This could be interpreted as indicating that the differing definitions of non-forest in k-NN and LCTS can explain half of the disagreements in that regard.Likewise, when the original k-NN quantitative rasters are used as reference, 63% of the pixels that are mixed forest in the SPOT map and that belong to a different forest type in the k-NN map would become mixed forest if k-NN used the same definition as SPOT.These two facts suggest that when two maps have different definitions for the same class, this will increase their level of disagreement in the allocation of that class; in the two cases above that we could quantify, it roughly doubled the number of pixels in disagreement.However, maps that have the same definition for non-forest, such as EOSD and k-NN, still show considerable disagreement: pixels with forest in one map and non-forest in the other add up to 177 Mha, which is even more than the 172 Mha of the same type of disagreement for the pair of maps (EOSD-LCTS) that have different definitions.Finally, uncertainties in the cross-walk of the original legends into our four classes, such as "young forest" in the original SPOT map or "treed wetland" in the original EOSD, can also contribute to additional disagreements between these maps and the rest, but their contribution should be minor given the way we dealt with them (see Section 2.3).
Discrepancies may also arise when the actual portion of terrain represented by each of the pixels occupying the same position in the rasters is not exactly the same.This would only affect the k-NN and LCTS maps, because the SPOT and EOSD maps were resampled to a much coarser pixel size and thus any misregistration at the native resolution will have a negligible effect on the 250 m pixels.Notwithstanding this fact, the resampling itself did cause a slight bias: it increased by 3% the proportion of needle-leaf forest in the SPOT map relative to the original 20 m pixel size, while the other classes decreased by roughly a percentage point.In the case of the EOSD, the resampling bias was even smaller: only mixed forest changed its proportion by more than 1 percentage point (+1.26).While resampling bias may certainly be larger locally, it cannot explain more than a small proportion of the disparities on the aggregate of millions of burned hectares (however, we acknowledge that other resampling methods such as nearest neighbor could have led to much different results, hence the effect of resampling should not be understated).For the two other maps that were not resampled, spatial mismatches and gridding artifacts [24] would still have no effect in the interior of a large homogeneous area covered by the same class, but in transitional or spatially heterogeneous areas, they certainly do.Spatial heterogeneity is a scale-dependent concept that in the context of land cover mapping can be defined as a configuration of the landscape fulfilling two conditions: (1) the actual patches, assuming they were accurately delineated using the same land cover classes but at a much finer spatial resolution, are similar in size to the pixels; and (2) in terms of composition, none of the classes predominate.In this situation, it is likely that adjacent pixels will display different classes, because the most abundant class in one pixel will seldom also be the most abundant in all of its neighbors.At the same time, slight offsets of the location of a pixel footprint could lead to a change of class.Therefore the k-NN and LCTS maps could disagree in spatially heterogeneous areas even in the absence of classification errors.
Past research has reported that spatial heterogeneity results in a greater number of mixed pixels that in turn decreases classification accuracy (e.g., [11,25]).Using sub-pixel fractional error matrices, Latifovic and Othlof [26] showed that the maximum achievable classification accuracy in spatially heterogeneous areas coincides with the average percent fraction of dominant land cover of pixels in those areas.From the analysis in the present study, we also found that spatial heterogeneity negatively affects the agreement among maps.The mean spatial heterogeneity (mean number of pixels in the eight-neighborhood with a different class, average of the four maps) of pixels where there is full agreement among all the maps was 1.2, whereas for pixels with no agreement it was 3.7 (for comparison, the mean value of the combined heterogeneity raster was 2.1, i.e., the average pixel has two neighbors from a different class).Conversely, the mean agreement for pixels in highly heterogeneous neighborhoods (more than six different neighbors) was 2.65 (lower than the mean agreement, 3.05), while for pixels surrounded by the same class it was 3.93 (i.e., full agreement for the great majority of those pixels).Therefore, as spatial heterogeneity increases, the level of agreement decreases.However, this relation is weak (R = ´0.53,p = 0.00).
Our conclusion from the previous paragraphs is that while differences in the definition of what is forest and what is mixed forest indeed contribute to some disagreements, and to a lesser degree so may the reclassification and resampling methods we applied, it seems that the majority of the discrepancies are likely due to classification errors.A cursory inspection of the confusion matrices of the six possible pairs of maps (populated using all pixels, Table 2) reveals that the overall agreement between different pairs of maps is similar, ranging from 59% (EOSD vs. k-NN and EOSD vs. LCTS) to 62% (EOSD vs. SPOT).This narrow interval suggests that all four maps must have a similar overall accuracy; otherwise, there would be one map for which the three confusion matrices that include the map would be much different than the other three that do not include it.At the class level, the variation in agreement is greater.The class with the best agreement is non-forest (which includes water), ranging from 49% (of non-forest pixels in EOSD that are also non-forest in k-NN) to 82% (the reverse in the same pair, i.e., 82% of non-forest pixels in k-NN are also non-forest in EOSD).The class with the worst agreement is mixed forest, ranging from 19% (of mixed forest pixels in LCTS that are also mixed forest in SPOT) to 57% (of mixed forest in SPOT that are also mixed forest in k-NN).However, even with a 57% agreement, the k-NN estimates of burned area belonging to mixed forest are, on average, five times greater than those of SPOT (Figure 3).When the confusion matrices were populated solely using pixels inside burned areas (not shown), the overall agreement between pairs of maps hardly changes (about two percentage points on average), decreasing in some cases (e.g., from 62% to 59% for SPOT vs. EOSD) and increasing in others (from 60% to 66% for k-NN vs. LCTS).Class-specific agreements change considerably more when only pixels in burned areas are considered.For example, the agreement in the mixed forest class for SPOT vs. EOSD goes from 35% to only 5%.This suggests that the commission and omission errors estimated from a confusion matrix from sample pixels across the entire map may not be good predictors of land cover mapping errors inside burned areas, because errors are not randomly distributed across the map and vary markedly from one place to another [22].Indeed, when we derived the confusion matrices from either the managed or the unmanaged forests, we obtained different overall agreements (the pair that changed the most was SPOT vs. EOSD: 59% in the managed forest compared to 65% in the unmanaged forest).Surprisingly, the agreement between forest classes was higher in the unmanaged forest than in the managed forest for all pairs of maps (12 percentage points better on average), even though the amount of reference data in the unmanaged forest is scant (e.g., k-NN had no training data there).The reason behind this apparent contradiction is that more than 2/3 of treed areas in the unmanaged forest are needle-leaf according to all maps, which increases the likelihood of pixel-wise agreement compared with the more diverse managed forest.Despite this, the largest unmanaged forest ecozone, the Taiga Shield, still shows large discrepancies among the maps in the breakdown of burned area by cover type (Figure 3).
One last point of interest was whether the observed overall agreement between the SPOT and LCTS maps (61%) was consistent with the overall accuracy (76% for the LCTS map and 79% for the SPOT map) that we derived for our four classes on the basis of the published confusion matrices of the two original classifications.The overall agreement between two maps is equal to the probability that a randomly chosen pixel has the same class in both maps.Assuming that errors in different maps are independent, this probability is the sum of the probability that both pixels are correctly classified (0.76 ˆ0.79) plus the probability that both pixels are incorrectly classified and end up with the same class ((1-0.76)ˆ(1 -0.79) ˆ0.36).Note that (1) 0.36 is the summation of the pairwise products of the proportion of each class in each map and (2) in reality, errors are positively correlated, thus the probability of simultaneous misclassification is greater than what we assumed, which makes our estimation conservative.The outcome of this calculation is 62%, almost equal to the actual figure of 61%.Since the other pair (EOSD and k-NN) has a similar agreement (59%), we can assume that their overall accuracy, if we had the original confusion matrices, would also be in the 75%-80% interval.Therefore, none of the maps seems better or worse than the others, and the most we can say, barring issues of comparability of the original legends and assessment methods, is that they all are relatively accurate compared with other large regional land cover mapping projects with reported accuracies ranging from 67% to 75% [18].

Conclusions
While all four maps agree that needle-leaf forest is the main component of burned areas in Canada, there is great disparity in the proportions of each broad land cover class (three forest types plus non-forest) that each map reports for the studied years.Proportions reported by one map can be more than double (or less than half of) those reported by another map, and the differences can be more than four times greater or smaller at the ecozone level.While differences among the maps in the definition of what is forest and what is mixed forest indeed contribute to some disagreements, the majority of the discrepancies seem to be due to classification errors, which vary from place to place and tend to be more frequent in spatially heterogeneous areas.Temporal or pixel-wise spatial mismatches between the maps seem to play a much smaller role.Given the observed discrepancies with just four classes and the lack of a clear candidate for the best map among the four datasets, it seems that none of the maps would be able to provide a sufficiently accurate answer to the question of what burns in Canada, at least one that would satisfy future reporting of GHG emissions for the unmanaged forest as envisioned in Canada's blueprint for forest carbon science [27].This is of concern, as more than half of the burned area occurs in the unmanaged forest (at least in 2011-2014), where the only full coverage datasets available are of the type we studied.Note that the managed forest relies on forest inventories of better quality than these data, although there are some data gaps.Therefore we concur with Kurz et al. [6] that further efforts are required to reduce data uncertainties in both the managed and unmanaged forest of Canada, of which improved land cover information is a necessary component (as a proxy for fuel type/loading, e.g., [28]).Better land cover information perhaps may be obtained at low cost by integrating existing datasets with new evidence, as was done by Tsendbazar et al. [29].In addition, there is a need to produce pixel-specific measures of accuracy, so that uncertainty can be propagated in spatial models that use land cover as a covariate [30].Tackling these deficiencies will help us to better understand how GHG emissions from Canada's vast unmanaged forests affect international commitments to curb global temperature increase.

Figure 1 .
Figure 1.The four land cover datasets in this study with the boundaries of the 10 ecozones.The pie charts represent the proportion of each class for each map across all the extent.

Figure 1 .
Figure 1.The four land cover datasets in this study with the boundaries of the 10 ecozones.The pie charts represent the proportion of each class for each map across all the extent.

Figure 2 .
Figure 2. Raster of agreement between the four maps.The pie chart represents the proportion of each agreement class (four classes, from no agreement to full agreement) across all the extent.Burned areas during 2011-2014 are overlain in yellow over the raster of agreement.The managed forest is the hashed area; the rest is unmanaged.

Figure 2 .
Figure 2. Raster of agreement between the four maps.The pie chart represents the proportion of each agreement class (four classes, from no agreement to full agreement) across all the extent.Burned areas during 2011-2014 are overlain in yellow over the raster of agreement.The managed forest is the hashed area; the rest is unmanaged.

Figure 3 .
Figure 3. Answers provided by the four maps to the question of what burned in Canada in each of the years 2011-2014.The amount of burned area in each year is indicated in parentheses below the year.Figure 3. Answers provided by the four maps to the question of what burned in Canada in each of the years 2011-2014.The amount of burned area in each year is indicated in parentheses below the year.

Figure 3 .
Figure 3. Answers provided by the four maps to the question of what burned in Canada in each of the years 2011-2014.The amount of burned area in each year is indicated in parentheses below the year.Figure 3. Answers provided by the four maps to the question of what burned in Canada in each of the years 2011-2014.The amount of burned area in each year is indicated in parentheses below the year.

Figure 4 .
Figure 4. Proportion of burned area in 2011-2014 in each ecozone according to each map.

Figure 4 .
Figure 4. Proportion of burned area in 2011-2014 in each ecozone according to each map.

Table 1 .
Definitions of the four land cover classes in each of the maps.

Table 2 .
Confusion matrices between the four maps based on all pixels.Figures inside the matrices are in millions of hectares (outside are percentages).Classes: NF, non-forest; NL, needle-leaf; BL, broadleaf; and M, mixed forest.