The central concept behind the dataset we present here is the combination, over a local moving window, of a static map of vegetation cover fractions with datasets of variables describing surface properties of vegetation that are retrieved from satellite observations. The methodological steps, including steps of masking, aggregating and cleaning, are summarized in the flowchart in Fig. 2. These steps are already described in our related work6, but the methodology is also reported here, along with some additional technical details, in order to facilitate the presentation of the dataset.

Figure 2: Schematic overview of the processing steps to generate the dataset. The different output files of the dataset correspond to the effects of vegetation cover change for: shortwave reflected radiation (SW ↑ ), longwave emitted radiation (LW ↑ ), latent heat (LE), sensible and ground heat fluxes combined (H+G), daytime and nighttime land surface temperature (LST d and LST n ), clear sky longwave emitted radiation (LW ↑ *) and albedo (α). Full size image

Data preparation

The input data for the methodology are of two types: surface property variables that change seasonally and static vegetation cover fraction maps. All have a common spatial resolution of 0.05. The seasonal variables are α, LE, LST d , LST n and LW ↑ *. They are all derived from measurements from the NASA Moderate Resolution Imaging Spectroradiometer (MODIS) instrument on-board the Aqua and Terra satellites. They are available at monthly temporal resolution, but a median value for each month is calculated from the years 2008 to 2012 to generate the 5-year climatology while retaining the seasonal cycle. The static vegetation cover fraction maps are derived from different satellite instruments and represent a single year: 2010. We describe each variable below, along with specific pre-processing each may have required.

Albedo

Albedo (α) is defined as the proportion of the incident light or radiation that is reflected by a surface. In this case, we are interested in the monthly average of the proportion of total radiation across the broadband shortwave spectrum reflected by the Earth's surface at 0.05 resolution. The NASA MCD43C3 albedo product provides 8-daily estimates of both directional hemispherical albedo (black-sky albedo) and bihemispherical albedo (white-sky albedo) based on multidate multispectral MODIS cloud-free observations collected over a 16-day moving window and a semi-empirical kernel-driven bidirectional reflectance model22. These white-sky and black-sky albedos correspond to theoretical situations in which incident radiation is either completely diffuse or completely direct. To obtain an estimate of real conditions without information on the fraction of diffuse radiation, we took the mean of both values. The mean is preferred here for its simplicity, since other techniques such as weighting or interpolation would require additional information to avoid any subjectivity. To obtain estimates at monthly temporal resolution, we selected those for which the 16-day period best corresponded with the 15th of each month. Data are available from the NASA LPDAAC website (https://lpdaac.usgs.gov/).

Latent heat flux

Latent heat flux (LE) is the flux of heat from the Earth's surface to the atmosphere that is associated with evaporation of water at the surface. The land component of this flux consists of evaporation of rain water intercepted by the canopy before it reaches the ground, transpiration through stomata on plant leaves and stems, evaporation from wet and moist soil and the sublimation of water vapour from ice and snow; and is otherwise known as terrestrial evapotranspiration (obtained by dividing latent heat flux by the latent heat of vaporization). The MOD16A2 product23 provides latent heat obtained by integrating several MODIS products (land cover, albedo, leaf area index, fAPAR) with meteorological data, delivered at 0.05 spatial resolution with monthly temporal resolution covering the regions from 60°S until 80°N. The data are available from the NTSG website (http://www.ntsg.umt.edu/project/modis/mod16.php).

Daytime and nighttime LST

These are the radiant temperatures of a surface measured respectively during the day and at night. The MODIS instrument on board the Aqua platform makes such measurements twice over its cycle at approximately 13:30 and 1:30 local time at the Equator. These times are close to those at which the minimum and maximum temperatures are expected. The MODIS land surface temperature algorithm provides such estimates at a monthly time step at 0.05 spatial resolution24,25 in the MYD11C3 product (we use collection 5 available from the NASA LPDAAC website https://lpdaac.usgs.gov/).

Surface upwelling longwave radiation

This is the outgoing infrared radiation emitted by the surface. A large part of this energy is absorbed by the atmosphere and later re-emitted towards the Earth (by clouds and greenhouse gases) or outwards to space. The upward longwave radiation (LW ↑ ) can be calculated from the surface temperature (T) and broadband emissivity (ε B ) using the Stefan-Boltzmann law:

(1) L W ↑ = ε B σ T 4

where σ is the Stefan-Boltzmann's constant (5.67×10−8 W m−2 K−4). We use LST d and LST n from the MYD11C3 product to estimate the mean surface temperature (T) over the entire day span using a simple average. The MYD11C3 product also provides emissivity estimates for various specific narrow bands in the middle and thermal infrared spectrum that can be used to obtain ε B using the empirical equation suggested by a dedicated study26:

(2) ε B = 0.2122 ε 29 + 0.3859 ε 31 + 0.4029 ε 32

where ε 29 , ε 31 and ε 32 are the estimated emissivities in MODIS bands 29 (8400-8700 nm), 31 (10780-11280 nm) and 32 (11770-12270 nm), respectively. Because the satellite can only measure during cloud-free conditions, it must be specified that the resulting monthly upwelling longwave radiation only refers to clear sky conditions, which we will denote using an asterisk: LW ↑ *.

Fractions of vegetation cover

These are derived from the 300 m global land cover map of the year 2010 provided by the European Space Agency's (ESA) Climate Change Initiative (CCI)27. A dedicated tool provided alongside the product enables users to transform the UNLCCS classification scheme as employed28 to continuous cover fractions at coarser scales via a ‘crosswalking’ table, typically to characterize the plant functional types used in land surface models29. We used this tool to generate two sets of maps for a set of 4 broad classes and another of 10 detailed classes (Table 1) based on the classes of the widely-used global vegetation classification scheme of the International Geosphere Biosphere Programme (IGBP). The crosswalking table we use to pass from UNLCCS to IGBP classifications is provided in Table 2 (available online only).

Table 1 Classes of vegetation cover used in the dataset. Full size table

Table 2 Cross-walking table used to translate land cover classes of the FAO LCCS system used in the ESA CCI maps to IGBP classes. Full size table

Retrieving the local biophysical signal of vegetation change

To identify the biophysical signal due to changes in vegetation cover we establish a relationship between vegetation cover fractions and the surface variables over a local moving window. As a result of this the direct biophysical effects of vegetation change considered here are local. This is valid both for the spatial extent of the cover change, which assumes at most a change of a complete fine resolution pixel (0.05°×0.05°), and for the origin of the change, i.e. we ignore indirect effects due to regional change from neighbouring regions. The moving window size is 5 by 5 pixel at 0.05° resolution, covering an area of approximately 25 km by 25 km over which the local climate is assumed to be uniform. To unmix the signal resulting from the mixed compositional land cover, for each window we apply a linear regression using a matrix X containing the vegetation fractions of each of the 25 pixels as explanatory variables and a vector y containing the 25 values of a given biophysical variable as response variable to obtain a vector of β coefficients:

(3) y = X β

This is equivalent to solving the following system of equations:

(4) { y 1 = β 1 x 11 + β 2 x 12 + .. . + β m x 1 m y 2 = β 1 x 21 + β 2 x 22 + .. . + β m x 2 m ⋮ y n = β 1 x n 1 + β 2 x n 2 + .. . + β m x n m

in which x ij represents the cover fraction of vegetation j in pixel i, for the n pixels in the moving window and the m classes that are considered. Once identified, we can use the β coefficients to predict the local y value corresponding to a given composition x, including that composed of a single vegetation cover j by setting x j =1 and all other x values to zero.

There is a problem, however, if the compositional predictor dataset X is used directly in the analysis. Compositional data can behave somewhat differently to ‘ordinary’, open or normal data, because compositions necessarily sum to one (for this reason they are also sometimes described as ‘closed’ data). Statistically, this can lead to spurious correlations between compositional components, and/or between compositional components and the response variable. Analysis of any given subset of compositional components can lead to very different patterns, results and conclusions30. Geometrically, all points defined by the compositions must fall in a simplex because their compositions sum to one. For a three part composition, this simplex is a triangular plane (i.e. it exists on a 2-dimensional surface). Whilst this compositional matrix has 3 columns, there are only (at most) 2 dimensions. A transformation of X is needed to reduce appropriately the dimensionality of this matrix for subsequent use in the regression.

The transformation we apply to reduce the dimensionality of X involves a singular value decomposition (SVD). This procedure is very close to a principal component analysis (PCA). The first step consists of centring all the columns of the predictor matrix X of vegetation fractions by removing the column means. We then apply the SVD:

(5) ( X − M ) = U D V t

where M is the appropriate matrix of column means, U and V are the matrices containing respectively the lefthand and righthand singular vectors, and D is a diagonal matrix containing the singular values (the standard deviations of the ensuing dimensions). Squared values of D indicate how much variance is explained by each (orthogonal) dimension. We implement a rule where as many dimensions from this SVD are retained as to conserve 100% of the original matrix's variation. In doing so, we reduce the dimensionality appropriately as described above, as well as remove what may be additionally redundant dimensions that can occur locally if, for instance, the only points in which 2 classes are represented have exactly the same values. To avoid having problems when there is too little or no information (e.g. if all pixels have exactly the same compositions), we added a pre-condition that there must be at least 10 pixels with different compositions. The final appropriately transformed predictor matrix of reduced dimension Z is then obtained by:

(6) Z = ( X − M ) V z

where the subscript z in V z indicates that the latter is composed of a subset of righthand singular vectors in V as selected from D as described above. The resulting predictor matrix Z can now be regressed onto the local biophysical variable y.

(7) y = Z β + ε

where Z has been augmented with a leading column of 1s to accommodate an intercept term in the regression. The standard manner to obtain an estimate of β is:

(8) β = ( Z t Z ) − 1 Z t y

Because the compositional predictor matrix X has been transformed to matrix Z, regression coefficients identified in the regression of Z onto y do not immediately provide information about the association between the various vegetation cover fractions and the surface property variables. In order to identify the z values associated with a particular vegetation (in that local analysis) we instead define a ‘dummy pixel’ whose composition contains only that vegetation class, with all other classes in the dummy pixel's composition set to zero. This pixel's composition is then transformed, and its y value predicted. This is the y associated with that vegetation type. Since we wish to do this for all compositional components of interest, we actually define a matrix P with as many rows as these compositional components that we wish to predict. P is centred on the same column means as above (M, specific to each local analysis), and then multiplied by the correct number of transposed right hand singular vectors (V z , again, specific to each local analysis).

(9) Z p = ( P − M ) V z

Predicted y p values for each vegetation type (identified by predicting the appropriately transformed ‘dummy pixels’) are then calculated as:

(10) y p = Z p β

The expected change in variable y associated with a transition from one vegetation type to another at the central pixel of the local window is then the difference between the y p predicted for each `pure' vegetation type:

(11) Δ y A → B = y B − y A

Beyond our primary interest in the change Δy for a given vegetation transition, we also assess the uncertainty associated with each of these differences. We consider uncertainty in terms of standard deviations, and thus, according to error propagation, the uncertainty for the difference due to the transition from A to B can be determined from:

(12) σ A → B = σ A 2 + σ B 2 − 2 σ A B

where σ A 2 and σ B 2 are the variances in the estimates of y for each vegetation type, and σ AB is their covariance. This covariance term is important as the uncertainties of the individually predicted z values are not independent given that they derive from the same regression model. The variances and covariances of all vegetation types can be obtained from the covariance matrix, which in turn is calculated as:

(13) Σ = Z p V a r [ β ] Z p t

The diagonal terms in Σ are the variances of individual predictions of (individual) vegetation classes. The off-diagonal parts of Σ hold the covariances between these predictions.

The whole procedure described above (variable transformation, regression and uncertainty estimation) is applied globally over 5 by 5 moving windows for the 3 biophysical variables for each of the 12 months of the year at 0.05° spatial resolution for each vegetation transition considered. The regressions are applied to data including information of all vegetation types (either 10 or 4, depending on whether the detailed or broad classification scheme is used) plus 4 non-vegetated classes (urban, water, snow/ice or bare soil), but predictions are only made for vegetated classes. Symmetric transitions yield identical results (e.g. Δy A→B =− Δy B→A ), and thus only a total of (102−10)/2=45 or (42−4)/2=6 transitions are calculated. The resulting maps only provide information for the pixels in which all 25 pixels in the moving window had information.

Masking out low vegetation co-occurrences

The method relies on there existing co-occurrences of vegetation classes within the local window. Furthermore, the statistical methods that are applied to these sets of points are more likely to provide reliable results when there are large and balanced presences of both vegetation classes of interest. We designed an indicator to quantify how two vectors describing the presences of two classes of vegetation, p A and p B , represent the potential range of compositional variation at the same time (i.e. how evenly balanced they are against each other, and how abundant they are with respect to the total composition). In a two dimensional space describing the presences of vegetation class A and vegetation class B, we consider a set of n points evenly distributed along the line B=1−A. These points represent the ideal situation and their positions are stored in vectors q A and q B . For each of these points q, we then calculate the Euclidean distance to all points in p, and retain only the smallest one that we store in a vector called d min :

(14) d m i n = min ( q A − p A ) 2 + ( q B − p B ) 2

The situation that would generate the largest cumulated distances is when all p points are located in the origin, i.e. when there is neither any of class A nor any of class B. This vector of maximum distance d max is defined as:

(15) d m a x = q A 2 + q B 2

The index of vegetation co-occurrence I c is defined as:

(16) I c = 1 − ∑ i d m i n , i ∑ i d m a x , i

I c ranges from 0 to 1, corresponding to a gradient of ‘no presence of either class’ to ‘full and evenly balanced presence of both classes’. For each pair of vegetation classes, a threshold of I c <0.5 is used to mask out from the results those pixels whose local windows do not provide enough co-occurrences. This threshold of I c =0.5 represents a situation in which the sum of the minimum distances between all points p and q is half of the distance between the origin and all q points.

Masking out high topographical variability

Restricting the analysis to a local window reduces the effect of major climatic gradients when comparing biophysical variables of neighbouring pixels. However, climatic gradients can also occur at much finer spatial scales as a result of vertical elevation change in a landscape, or topographical relief. To factor out such a potentially confounding effect, we created a mask to remove areas in which neighbouring pixels should not be compared because of within pixel differences in relief. This mask is constructed based on fine-spatial resolution digital elevation models: the gap-filled SRTM 90 m Digital Elevation Database v4.131 for the region it covers (between 60°S and 60°N) and GTOPO30 at 1 km for the rest of the globe (data available from the U.S. Geological Survey and distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at USGS/EROS, Sioux Falls, SD. http://lpdaac.usgs.gov). The mean and standard deviation of the elevation is calculated for each 0.05° by 0.05° grid cell to produce two data layers: μ h and σ h . These two layers are used to derive three indicators of local topographical relief using the values of all i pixels in the 5 by 5 pixel moving window. The first, v 1 , is just the average standard deviation across the the moving window:

(17) v 1 = 1 n ∑ i = 1 n Σ h , i

High values of v 1 indicate hilly terrain over the considered scale, which should be discarded from the analysis. The second index, v 2 , indicates how different the mean elevation within the central pixel is from the average elevation in the local window:

(18) v 2 = | µ h − 1 n ∑ i = 1 n µ h , i |

High values of v 2 can further isolate unwanted pixels that may not be identified using v 1 alone. The third, v 3 , is similar to v 2 but compares the central pixel's standard deviation of elevation with the standard deviation across the moving window:

(19) v 3 = | σ h − v 1 |

High values of v 3 can isolate the odd undesirable pixel whose within-pixel elevation has high variance but an average value close to the average height across the local window. These three indicators are combined together in a single layer depicting all pixels satisfying all of the following conditions: v 1 <50 m, v 2 <100 m and v 3 <100 m. Pixels that fail any of these conditions are masked out from all the layers of results.

Spatial aggregation

The maps resulting from the local space for time analysis need to be spatially aggregated from 0.05° to 1° grid cells to be used alongside data from the CERES instrument, which provides the information necessary to close the energy balance. Aggregating to 1° also has other advantages, namely: (1) a mean difference of a variable associated with change from one vegetation type to another may be assumed to be more accurate than any individual estimate at finer scale; (2) this scale is simpler to map and visualize at global level; and (3) it is more comparable to results from land surface models. Because each 0.05° estimate of Δy includes an associated estimation of its uncertainty, this uncertainty can be used to down-weight less reliable values during the aggregation procedure. The typical approach to do so is weighting based on the inverse of the uncertainty:

(20) Δ y - = ∑ i Δ y i / σ i 2 ∑ i 1 / σ i 2

where Δy̅ is the mean aggregated value, whose uncertainty is calculated as:

(21) σ Δ y - 2 = 1 ∑ i 1 / σ i 2

However, these formulations do not account for the spatial auto-correlation generated by the moving window (1 to 20 pixels may be common between two nearby estimates depending on the possible overlap of their respective 5 by 5 windows). This auto-correlation problem may be compounded further when only a clustered set of 0.05° samples are available within the 1° by 1° area. This can occur due to the topographical masking, or because two vegetation types only co-occur over a small part of the 1° grid cell.

To tackle this auto-correlation, we employ a more generic weighting approach. The weights depend not only on the uncertainties estimated from the regressions as above, but also on how each window is correlated with every other window within the area of 1° grid cell. This information is summarized in a 400 by 400 matrix R a containing the fraction of overlap between every pair of windows. The information in R a is combined with that of the pixel-wise uncertainties that are embedded in D a , a diagonal matrix containing the uncertainties in its diagonal, to build a covariance matrix Σ a (the subscript a is used to differentiate these matrices involved in this aggregation step from those used before):

(22) Σ a = D a R a D a t

The vector of weights is then obtained as:

(23) w = 1 1 t Σ a − 1 1 Σ a − 1 1

which can then be used to calculate the aggregated Δȳ as:

(24) Δ y ̄ = ∑ i w i Δ y i

while the aggregated uncertainty σ Δȳ 2 is given by:

(25) σ Δ y ̄ 2 = w t Σ a w = 1 1 t Σ a − 1 1

When the windows have no auto-correlations, both equations 24 and 25 fall back to the simpler weighting formulas of equations 20 and 21. The aggregation procedure is applied to all data layers.

Detection and treatment of outliers

Despite all efforts to characterise uncertainty and reach representative values, the results can still contain unrealistic values. The main reasons for this might be that uncertainties in the input data (the remote sensing surface property variables and the vegetation cover fraction maps) are not explicitly taken into account. As a final step to remove possible outliers, we remove all values for grid cells in which there are not at least 20 samples at 0.05° spatial resolution. Lastly, we also remove values that are statistical outliers based on the distribution of the entire dataset. All data layers are available with their associated uncertainty.

Closing the surface energy balance

The local unmixing step can only be applied to those variables available at the 0.05° spatial resolution (namely α, LE, LW ↑ *, LST day and LST night ), meaning some components of the surface energy balance are missing. The full surface energy balance is expressed as:

(26) S W ↓ − S W ↑ + L W ↓ − L W ↑ = H + L E + G

SW ↓ , SW ↑ , LW ↓ and LW ↑ are respectively the downwelling and upwelling radiative fluxes in the shortwave or longwave parts of the spectrum, LE is the latent heat flux, H is the sensible heat flux and G is the ground heat flux. We derive the terms of the energy balance combining MODIS-based datasets with the EBAF-Surface Product derived from the NASA Clouds and the Earth's Radiant Energy System (CERES) instrument. This dataset (CERES EBAF-Surface Ed2.8) provides a closed and gap-filled surface energy balance at 1° that is consistent with CERES top-of-atmosphere irradiance measurements32. For the specific goals of this analysis we are interested in how the terms of this equation change according to a change in vegetation cover, i.e.:

(27) Δ S W ↓ − Δ S W ↑ + Δ L W ↓ − Δ L W ↑ = Δ H + Δ L E + Δ G

We make the assumption that the changes in vegetation cover that are considered here are too small (i.e. maximum 0.05°) to generate strong feedbacks in the cloud regime, and as a consequence we assume ΔSW ↓ =0 and ΔLW ↓ =0. The change in reflected shortwave radiation can be expressed in terms of albedo (α) and incoming shortwave radiation (ΔSW ↑ =Δα×SW ↓ ), the latter being available from CERES data at 1°. Although we derived estimates of changes in upwelling longwave flux satellite measurements at 0.05°, these refer to clear-sky conditions only (i.e. when the satellite instrument can measure the ground unobstructed by clouds) while other fluxes are representative of all cloud conditions. As a proxy for the effect of cloudiness, we used a correction factor based on the ratio of all sky (LW C↑ ) to clear sky (LW C↑ *) longwave upwelling estimated by CERES (ΔLW ↑ = (ΔLW C↑ /LW C↑ *)×ΔLW ↑ *, where the asterisk indicates values for clear sky conditions). By re-writing and simplifying the equation above, the expression describing the change in the residual flux, composed of both sensible and ground heat fluxes, becomes:

(28) Δ ( H + G ) = − ( Δ α ) S W ↓ − ( Δ L W C ↑ / L W C ↑ ) × Δ L W ↑ − Δ L E

We apply this expression to every 1° pixel for every month of the time-series and every vegetation transition based on the previously calculated datasets of Δα, ΔLW ↑ * and ΔLE. To have all terms of the energy balance on equal footing and with the same sign convention, we also explicitly produced datasets of shortwave reflected radiation (ΔSW ↑ ) and full-sky longwave emitted radiation (ΔLW ↑ )

Code availability

Most of the processing has been done using R version 3.3.2. The code can be made available upon request on a case by case basis. The aggregation and transformation of the ESA CCI land cover maps is done using the dedicated user tool (version 3.12) that can be accessed at: http://maps.elie.ucl.ac.be/CCI/viewer/download.php#usertool

Data Records

The entire dataset (Data Citation 1) is composed of 16 separate netCDF files, as listed in Table 3 (8 variables, for each of two vegetation cover classifications). The 8 variables are: (1) those based only on MODIS data (LE, albedo, LWsfc, LSTday and LSTnight); (2) those requiring a combination of one of the previous along with CERES data (SWreflected and LWemitted); and (3) those requiring a combination of all the former (HG). The vegetation classifications are either detailed (IGBPdet, 10 classes) or generic (IGBPgen, 4 classes) as described in the methods section and in table 1. The names of the netCDF files include both the variable and the land cover classification scheme (e.g. HG_IGBPdet.nc for Δ(H+G) given the detailed land cover scheme, and that therefore includes 45 transitions).

Table 3 List of separate files contained in the dataset. Full size table

All files provide information across the global domain at a spatial resolution of 1°×1° and monthly temporal resolution over a climatological year representing the period 2008-2012. Each netCDF file has the same structure based on 4 different dimensions: latitude (lat), longitude (lon), month (mon) and vegetation transition code (iTr). Three dependent variables are recorded over these 4 dimensions for each variable Z: the difference in variable Z for a given vegetation transition (Delta_Z), the uncertainty (provided as a standard deviation) associated with the estimation of this difference in Z (SD_Delta_Z), and the number of fine spatial resolution (0.05°×0.05°) samples from which each aggregated (1°×1°) value is derived (N_Z).

The vegetation transition code (iTr) is a two or three digit number, representing the class of origin (first digit), and the class of destination (second, or second and third digit), with classes coded as in table 1. For example, for the 10 class system IGBPdet, the code iTr=19 represents the transition from class 1 (EBF) to class 9 (CRO), while the code iTr=23 represents the transition from class 2 (DBF) to class 3 (ENF). Because there are 10 classes, any transition from classes 1-9 to class 10 will result in a three digit transition code (e.g. iTr=110 represents the transition from class 1 (EBF) to class 10 (wetlands, WET)). There is no ambiguity with a three digit transition code, they are always from a (1 digit) class to class 10. Transitions from a larger number to a smaller number are not encoded as they are equal to the transition from the smaller number to the greater number, albeit with a change of sign (e.g. iTr=32 would be equivalent to iTr=23 with a negative sign added to the mean change).

Technical Validation

The dataset successfully reproduces the major patterns that are expected following vegetation cover change. Figure 3 reports the mean annual effect of deforestation for each of the 8 variables in the dataset, summarized across climatic gradients of mean annual temperature and annually cumulated precipitation (both from CRU data v4.00 at 0.5°× 0.5° resolution). The change in albedo (and shortwave reflected radiation) is a consistent brightening in all climates (since forest are generally darker than crops or grasses1), but with much higher values in cold climates, where the snow effect is expected to exacerbate the change, and very dry areas, where grasses remain dry for large parts of the year or crops are not so densely sown and crop seasons are shorter, allowing a stronger influence of the bright dry soil on the signal. The change in longwave emitted radiation and daytime LST is stronger in warmer climates where an increase is expected caused by the strong decrease in latent heat, which in turn is caused by the fact that herbaceous plants typically have shallower roots and less access to water than trees. The changes in nighttime LST are milder than those during daytime, but they do show a change in sign, with a decrease in colder climates and a slight increase in warmer ones. The reason behind this is probably the reduction in surface roughness going from trees (particularly needleleaf which dominate in colder climates) to crops/grasses, causing a stable stratification of the air at night, which is particularly evident in colder climates and at lower levels of radiation. The changes in the residual fluxes (H+G) are the consequence of all the other fluxes in the surface energy balance.

Figure 3: Summary of annual changes in all variables due to a given type of vegetation cover change plotted in a climatic space. Values of a given biophysical variable have been divided by the standard deviation of the entire population prior to aggregation in climate space. The climate values come from CRUNCEP and cover the period 2008–2012. Full size image

A validation with ground-based measurements of surface energy balance fluxes would be desirable to make a proper detailed evaluation of the present dataset. However, no such network of measurements currently exists. The closest candidate are the flux-tower measurements from FLUXNET2015, but there are insufficient sites with comparable climate conditions and contrasting vegetation types to cover the spatial and thematic (i.e. vegetation type transitions) extents in the datasets. In the absence of proper validation, an indicative evaluation is proposed based on a selection of paired FLUXNET sites that satisfy similarity criteria based on ERA-Interim reanalysis data at each site. These criteria require that to be paired, sites must have: (1) less than 10% difference in cumulated annual precipitation; (2) less than 2°C difference in mean annual air temperature; and (3) similar degree of continentality, defined as less than 0.5° difference in standard deviation of monthly air temperature. Only paired sites including transitions between forests and grasses or crops are considered, resulting in 9 pairs whose difference in turbulent fluxes are compared with the values from the present dataset at each site in Fig. 4 (radiative fluxes were not sufficiently available in the dataset). These comparisons are not expected to match closely as the spatial supports are radically different: the flux-tower measurements represent a specific patch of area of less than 1 km2 while the remote sensing data pixels cover the average conditions over 1°×1°. However, the seasonal profiles of our datasets do generally show a reasonable degree of similarity with those from the paired sites, and this similarity increases when they are aggregated according to general classes as shown Fig. 5.

Figure 4 Comparison of the dataset against individual paired sites of flux-tower measurements. Full size image

Figure 5: Comparison of the dataset against averaged paired sites of flux-tower measurements. The value of N indicates the number of paired sites involved in the average. Full size image

Usage Notes

Several assumptions were necessary when constructing this dataset, and should be considered when interpreting or reusing it. An important point to emphasise is that, although the data layers are provided with a spatial resolution of 1°, the values themselves represent a spatial average of effects that would occur for vegetation cover change at the scale of 0.05° or less, and thus represent direct local effects. A total change of vegetation cover over an area of 1°×1°, which is highly unrealisitic, would potentially cause a much larger perturbation involving both direct and indirect effects.

We also assume vegetation cover is the only driver of changes in biophysics within the local moving window of 0.25°. Areas with strong elevation gradients are masked out to filter topographic effects, but other changes, in soil properties for example, are disregarded. To close locally the surface energy balance, we also need to assume that changes of vegetation cover within the 0.05°×0.05° area do not generate indirect cloud feedbacks with the atmosphere that would further change the energy balance. For example, due to their higher roughness forested landscapes can induce more cloud formation than grasslands33, and these small clouds can change the energy balance by reducing incoming shortwave radiation, increasing incoming longwave radiation and exporting latent heat. Differences in soil moisture under contrasting vegetation covers could also contribute to indirect changes in the surface energy balance that cannot be resolved in this dataset. Assessing these kind of effects would require the use of land surface models coupled with dynamic regional climate models which are beyond the scope of this work. The assumption to neglect these feedbacks relies on the fine scale of the analysis and on the typical lateral movement of air masses due to wind that ultimately advect air masses to different grid cells.

The dataset is based on satellite observations obtained during the period 2008–2012. These values should remain valid for various years in the past or the future if the background climate does not change substantially9. Changes in background climate, such as an average increase in temperature that reduces the length of the snow cover period or intensifies heat stress in summer, will change the albedo and the evapotranspiration signal respectively. User are not encouraged to project these values in scenarios of future climates without considering carefully the associated changes in climate at every pixel.