Automatic Generation of Land Surface Emissivity Maps

1. Introduction

The remote sensing measurement of the land surface temperature (LST) from satellites provides an overview of this magnitude on a continuous and regular basis. The study of its evolution in time and space is a critical factor in many scientific fields such as weather forecasting, detection of forest fires, climate change, etc.

The main problem of making this measurement from satellite data is the need to correct the effects of the atmosphere and the land surface emissivity (LSE). Nowadays, these corrections are usually made using a split-window algorithm, which has an explicit dependence on land surface emissivity.

Therefore, the aim of our work was to define an enhanced vegetation cover method and develop a computer system that used it, in order to calculate and generate, automatically, maps of land surface emissivity from images of the AATSR (Advanced Along Track Scanning Radiometer) onboard the ENVISAT satellite.

The most innovative part of our method is that we provide it with the resources and the capability to calculate the most accurate coefficients according to the specific characteristics of each area (vegetation cover fraction, vegetation type, season, etc.). This allows the method to be applied to generate large-scale maps of this magnitude (Caselles et al., 2009).

On the other hand, the current procedure (global, fully operational and supported by ESA) for obtaining the emissivity from an AATSR pixel (Noyes et al. 2007) causes systematic errors when calculating the temperature of 2 to 5 K (Coll et al. 2005), showing that the current classification and the vegetation cover maps made with a resolution of 0.5º x 0.5º could be highly improved and provided with the same spatial resolution of the AATSR images (1km x 1km). This is the main reason of this paper.

In this chapter, this new method is presented, with its algorithm, and it is applied to several different types of vegetation in AATSR images of Europe, making all the calculations automatically with the developed software.

Eventually, an on field validation of the method was carried out by comparing the data of the generated emissivity maps (as the one in figure 6) with the values obtained in previous campaigns (Coll et al. 2005) carried out in the area of rice fields of Valencia, Spain (Caselles et al., 2009). An error of less than ±0.01 in the land surface emissivity assessment was successfully obtained.

Its validation was made by comparing the obtained results and the values measured in previous field campaigns carried out in the area of rice fields of Valencia, Spain.

2. Methodology

So far, different methods to obtain the land surface emissivity based on different ideas have been purposed (Caselles et al., 1997). The main disadvantages they present, especially when your final aim is to calculate this magnitude automatically, are:

1. Their high complexity, since it is difficult to apply them operatively, carrying out huge computational calculations.

2. Their error, because the propagation of errors is bigger in complex algorithms.

3. Their bias, which can be introduced if the approaches of the model are not fulfilled exactly.

For that reason, an enhanced mathematically simple method (Caselles et al., 2009) was defined, without important bias, to obtain the surface emissivity from satellite images, inspired by the results of Valor & Caselles (1996), who purposed a relation between thermal emissivity and the Normalized Difference Vegetation Index (NDVI).

The general philosophy of this new method is similar to some extent to the algorithm for LSE estimation used by the LST product of Terra-MODIS (Snyder et al. 1998), and for LSE estimations in Meteosat-SEVIRI (Peres & DaCamara 2005; Trigo et al. 2008)

2.1. Enhanced land surface emissivity model

As we have already mentioned, to produce the emissivity maps, we used an improved geometric model based on the one described in Valor & Caselles (1996). This model can be summarized in the following equation (1).

$$\epsilon ={\epsilon }_v{P}_v+{\epsilon }_g(1-P_v)+4<d\epsilon >P_v(1-P_v)$$E1

where ε_v and ε_g are, respectively, the vegetation and soil emissivities, <dε> is the effective cavity term and P_v is the vegetation cover fraction. It allows us to calculate the effective emissivity in a heterogeneous surface from a land use map and vegetation cover fraction image.

Therefore, we need to know the pure surfaces emissivities (ε_v and ε_g), that is, the emissivity of the vegetation and the existing ground under it, as well as the effective cavity term (<dε>). Since we do not have field measures from all over Europe and our system is aimed to be applied to the whole planet, the only existing possibility is to use average values obtained considering the possible variation ranges, experimentally observed.

In order to estimate the vegetation cover fraction (P_v), if we know the reflectivity values and i_v and i_g are the NDVI values obtained for a full vegetated surface and for a bare soil one, respectively, and K is given by:

$$P_v=\frac{1-\raisebox{1ex}{$i$}\!\left/ \!\raisebox{-1ex}{$i_g$}\right.}{\left(1-\raisebox{1ex}{$i$}\!\left/ \!\raisebox{-1ex}{$i_g$}\right.\right)-K\left(1-\raisebox{1ex}{$i$}\!\left/ \!\raisebox{-1ex}{$i_v$}\right.\right)} K=\frac{{\rho }_{2v}-{\rho }_{1v}}{{\rho }_{2g}-{\rho }_{1g}}$$E2

being ρ_2v and ρ_1v the near infrared and red vegetation reflectivities, respectively, and ρ_2g and ρ_1g the same measurements made on bare soil.

2.2. The land cover classification

In order to be able to obtain the most accurate values for the coefficients that have some dependence on vegetation and ground properties, we made our own collection, because we needed a concrete and adapted classification to the purpose of this chapter.

For that purpose, we studied some of the most famous existing land cover classifications, specially the Corine Land Cover (Buttner et al. 2004) of the European Environment Agency and the Ionia GLOBCOVER (Bicheron et al. 2008) of the European Space Agency. It is possible to read a complete analysis of both of them in Neumann et al. (2007). Eventually, we decided to use the second one, obtained from MERIS images, as our starting point, since it was updated recently to a newer version 2.2 (year 2008) and its higher spatial resolution (300 meters).

Furthermore, we made our own collection (see Table 1), by grouping the GLOBCOVER classes with similar emissivity characteristics (see Table 2), because we needed a more concrete and adapted to the purpose of this work classification. We did this following the land cover classification methodology explained by the FAO in DiGregorio & Jansen (2000).

Emissivity Class	AATSR-11 μm	AATSR-12 μm
Flooded vegetation/ crops/grasslands	εv =0.983±0.005	εv =0.989±0.005
	εg =0.970±0.005 (ground)	εg =0.977±0.004 (ground)
	εg =0.991±0.001 (water)	εg =0.985±0.001 (water)
	<dε"/>=0	<dε"/>=0
Flooded forest/shrubland	εv =0.981±0.008	εv =0.982±0.009
	εg =0.970±0.005 (ground)	εg =0.977±0.004 (ground)
	εg =0.991±0.001 (water)	εg =0.985±0.001 (water)
	<dε"/>=0.014±0.004 (ground)	<dε"/>=0.010±0.003 (ground)
	<dε"/>=0.004±0.001 (water)	<dε"/>=0.007±0.002 (water)
Croplands/grasslands	εv =0.983±0.005	εv =0.989±0.005
	εg =0.970±0.005 (ground)	εg =0.977±0.004 (ground)
	<dε"/>=0	<dε"/>=0
Shrublands	εv =0.981±0.008	εv =0.982±0.009
	εg =0.970±0.005 (ground)	εg =0.977±0.004 (ground)
	<dε"/>=0.014±0.004 (ground)	<dε"/>=0.010±0.003 (ground)
Broadleaved/needleleaved deciduous	εv =0.973±0.005	ev=0.973±0.005
forest	εg =0.970±0.005 (ground)	εs =0.977±0.004 (ground)
	<dε"/>=0.019±0.006	<dε"/>=0.015±0.004
Broadleaved/needleleaved evergreen	εv =0.989±0.005	εv =0.991±0.005
forest	εg =0.970±0.005 (ground)	εg =0.977±0.004 (ground)
	<dε"/>=0.019±0.005	<dε"/>=0.015±0.004
Urban area	ε =0.969±0.006	ε =0.976±0.004
Bare rock	ε =0.93±0.05	ε =0.95±0.05
Water	ε =0.991±0.001 (water)	ε =0.985±0.001 (water)
Snow and ice	ε =0.990±0.004	ε =0.971±0,014

Table 1.

Emissivity classes with the values for the parameters of the method in the 11 and 12 μm channels.

Emissivity Class	GLC Class	GLC Label
Flooded	11	Post-flooding or irrigated croplands (or aquatic)
vegetation/crops/grasslands	13	Post-flooding or irrigated herbaceous crops
	180	Closed to open ("/>15%) grassland or woody vegetation on regularly flooded or waterlogged soil - Fresh, brackish or saline water
	185	Closed to open ("/>15%) grassland on regularly flooded or waterlogged soil - Fresh or brackish water
Flooded forest/shrubland	170	Closed ("/>40%) broadleaved forest or shrubland permanently flooded - Saline or brackish water
Croplands/grasslands	14	Rainfed croplands
	15	Rainfed herbaceous crops
	20	Mosaic cropland (50-70%) / vegetation (grassland/shrubland/forest) (20-50%)
	21	Mosaic cropland (50-70%) / grassland or shrubland (20-50%)
	120	Mosaic grassland (50-70%) / forest or shrubland (20-50%)
	140	Closed to open ("/>15%) herbaceous vegetation (grassland, savannas or lichens/mosses)
	141	Closed ("/>40%) grassland
	150	Sparse (<15%) vegetation
	151	Sparse (<15%) grassland
Shrublands	16	Rainfed shrub or tree crops (cash crops, vineyards, olive tree, orchards…)
	30	Mosaic vegetation (grassland/shrubland/forest) (50-70%) / cropland (20-50%)
	130	Closed to open ("/>15%) (broadleaved or needleleaved, evergreen or deciduous) shrubland (<5m)
	131	Closed to open ("/>15%) broadleaved or needleleaved evergreen shrubland (<5m)
	134	Closed to open ("/>15%) broadleaved deciduous shrubland (<5m)
	152	Sparse (<15%) shrubland
Broadleaved/needleleaved deciduous forest	40	Closed to open ("/>15%) broadleaved evergreen or semi-deciduous forest ("/>5m)
	41	Closed ("/>40%) broadleaved deciduous forest ("/>5m)
	50	Closed ("/>40%) broadleaved deciduous forest ("/>5m)
	60	Open (15-40%) broadleaved deciduous forest/woodland ("/>5m)
	90	Open (15-40%) needleleaved deciduous or evergreen forest ("/>5m)
	91	Open (15-40%) needleleaved deciduous forest ("/>5m)
Broadleaved/needleleaved	32	Mosaic forest (50-70%) / cropland (20-50%)
evergreen forest	70	Closed ("/>40%) needleleaved evergreen forest ("/>5m)
	92	Open (15-40%) needleleaved evergreen forest ("/>5m)
	100	Closed to open ("/>15%) mixed broadleaved and needleleaved forest ("/>5m)
	101	Closed ("/>40%) mixed broadleaved and needleleaved forest ("/>5m)
	110	Mosaic forest or shrubland (50-70%) / grassland (20-50%)
Urban area	190	Artificial surfaces and associated areas (Urban areas "/>50%)
Bare rock	200	Bare areas
	201	Consolidated bare areas (hardpans, gravels, bare rock, stones, boulders)
	202	Non-consolidated bare areas (sandy desert)
	203	Salt hardpans
Water	210	Water bodies
Snow and ice	220	Permanent snow and ice
	230	No data (burnt areas, clouds,…)

Table 2.

Correspondence between GLOBCOVER classes and the emissivity classes.

A reduction from the 22 initial classes of the GLOBCOVER classification (Bicheron et al. 2008) was carried out to 10 classes, taking into account similarities between related classes from the point of view of its components and their typical structure. So, for each vegetated area we calculated average values of the ground (ε_g) and vegetation emissivity (ε_v) in each spectral band (11 and 12 μm) from the spectra of soils and vegetation emissivity given in the ASTER spectral library version 2 (Baldridge et al., 2009).

Along with these coefficients, it has been calculated an average value of the cavity term (<dε>) too, taking into account the structure of each vegetation type as described in the GLOBCOVER classification (Bicheron et al., 2008), by using the procedure defined in Valor & Caselles (2005).

In the case of water surfaces, snow and ice, or bare soil it has been allocated directly emissivity values from samples of the ASTER library. On the other hand, for urban areas, it has been used the effective value proposed by Valor et al. (2000), determined by the emissivity values of urban materials (mainly concrete, asphalt, ceramics) and the structure of buildings. Table 1 shows the defined classes, their descriptions, and the values applicable to the equation (1).

2.3. The automatic developed system

The software developed to produce the emissivity maps uses the land cover map, a table with the information for each land cover class and the different AATSR images of the area we want to produce the map. So, it extracts the required data from them, in order to be able to apply the mathematical model (1), previously explained.

These AATSR images have its own format, defined in the specification of the sensor. The system obtains the coordinates of the studied surface and the measured values at each channel of the sensor.

The land cover classification map used is similar to the GLOBCOVER v2.2 one (see Bicheron et al. (2008)). It uses the GeoTIFF format (as detailed in Richter & Ruth (2000)) and the system uses it to obtain the ground type of a pixel given by the coordinates of an AATSR image. This image has a spatial resolution of 300 m, while the ones produced by the sensor have a 1 km resolution. That is the reason why the system has been equipped with an interpolation by proportion of occupied areas algorithm to be able to combine them accurately, knowing the geographical coordinates of both images or maps.

A detailed example of the interpolation that this algorithm does between an AATSR pixel (red square) and its correspondent GLOBCOVER pixels (black edge squares) is shown in figure 1. Since the algorithm know the geographical coordinates of the center of the AATSR pixel and the geographical coordinates of the GLOBCOVER pixels, it knows where the surface of this pixel is in the GLOBCOVER image. In the example, it is centered in the coordinates (15.75, 20.32). Then, knowing the resolution of both pixels, the algorithm will detect which GLOBCOVER pixels are completely or partially inside by the AATSR one and the actual proportion of each of the firsts is occupied by the second one.

Figure 1.

Example of an AATSR pixel (1 km) in red interpolated with GLOBCOVER pixels (300 m) with black edge, according to the interpolation by proportion of occupied areas algorithm. Numbers between parentheses represent coordinates where the AATSR pixel is located in the GLOBCOVER image. Numbers near the green lines are the portion of the size of each side of every GLOBCOVER pixel inside the AATSR one.

The proposed interpolation algorithm obtains the different types of vegetation and soil that form each AATSR pixel and very accurately estimates the proportion of the area that each vegetation and soil type represents. First, the pixels of 300 m that the 1 km one overlaps with are obtained, identifying the GLC class each of them belongs and therefore the emissivity coefficients associated with them.

Subsequently, the exact area of each GLC pixel overlapped by the AATSR pixel is calculated based on their geographical coordinates and resolutions. eventually, the algorithm is able to estimate the values that need to be applied to each coefficient of equation (1), calculating each one as the weighted average of the values for that coefficient related to all the emissivity classes involved, and determining the influence of each emissivity class by the percentage that its area represents in the total area occupied by the AATSR pixel.

Therefore, the system processes all the pixels for each AATSR image, one by one, using the following algorithm (see flowchart in figure 2), reading the reflectivities (in the red and infrared channels) of each pixel and applying to them the mentioned model (1) to obtain the emissivity. Once all the pixels of one AATSR image are processed, a map with the calculated emissivities is generated for the same original surface studied by the sensor.

Finally, the system produces an output file, also following the mentioned GeoTIFF format, where it is stored: the average emissivity map, a confidence band, a land cover map, one NDVI map and one vegetation cover fraction (P_v) map, all of them for the original AATSR studied area. Any interested reader may ask the authors for a copy of this program, if desired.

Figure 2.

a) Main flowchart of the system designed to produce emissivity maps. b) Main flowchart of the system designed to produce emissivity maps.

3. Results

Once the enhanced method to calculate the land surface emissivity and the algorithm used by the system developed are defined. In this chapter, they are applied to Europe, in order to obtain a map of mean emissivities to a particular month of summer. Subsequently, a validation of that image is carried out to evaluate the system, comparing its results with actual measures.

3.1. Emissivity maps for Europe

As a result of using the above explained system, we have produced the emissivity map for Europe shown in Figure 3. It has been generated combining the output products of the system (in GeoTIFF format) for a set of 183 AATSR images, measured by this sensor in July, 2007. So, for each AATSR image, we have obtained the corresponding emissivity product and finally, we have joined them to create this composite. Also, figures 4 and 5 give the vegetation cover fraction and the confidence band for the same European area, respectively.

Figure 3.

Average emissivity map (between channels 11 and 12 μm) for Europe (July, 2007).

Figure 4.

Vegetation cover fraction (P_v) map for Europe (July, 2007).

Figure 5.

Confidence band for the product (valid pixels are in black).

3.2. Validation with field measurements

The validation of the whole system was carried out by comparing the data of the generated emissivity maps (as the one in figure 6) with the values obtained in previous campaigns (Coll et al. 2005) carried out in the area of rice fields of Valencia, Spain (Caselles et al. 2009). In order to compare them, we took all the measures data and the corresponding (by geographical coordinates) emissivity pixels and calculated the mean and the range of error for all, in each channel, as shown in Table 3.

Figure 6.

Average emissivity map (between channels 11 and 12 μm) for the AATSR sensor of the validation area (Valencian Community, Spain, 20/07/2007).

Channel (μm)	Measured value	System’s value
10.5-11.5	0.985 ± 0.002	0.982 ± 0.001
11.5-12.5	0.980 ± 0.005	0.988 ± 0.002

Table 3.

Comparison between experimentally measured values of emissivity and the ones obtained by the system for the same area of rice fields (Valencia, Spain).

Comparing both values, a difference of less than 1% (see Table 3) is obtained. This difference is derived from the original model, as explained in Valor & Caselles (1996), and it represents an error of ±0.5 K (Mira et al. 2007) when studying the temperature.

Since the magnitude of the error of this model is very important because the emissivity error determines the temperature error and the temperature is a determinant input parameter in a wide number of models, such as circulation models, global change models, energy balance models, climatic models, etc. Therefore, our present research objective is to work to keep reducing the magnitude of this error in the future. Some possible lines to investigate are the emissivity angular dependence and the emissivity soil water content dependence.

4. Conclusions

In this work, we have developed a system able to obtain, automatically, the surface emissivity with an ±1% error (as explained in Valor & Caselles (1996)). Although we have validated it in the area of Valencia, so we have verified that commits this error; in the future we hope to be able to validate it with field measures of other European regions.

Our algorithm combines, automatically, the data measured by the sensor onboard a satellite with the information contained in the land cover classification maps. We believe this algorithm, as the developed system that implements it, come to fill the gap that nowadays exists between the different methods and techniques to generate emissivity maps from satellite images.

This new model, by using the simple interpolation function between areas measured by the satellite and the land cover classification maps, could be the best solution, since it is easy to understand that finding the correct type of soil and vegetation for each area, in order to obtain and apply the most appropriate coefficients for it is always a more accurate procedure than considering a single value of emissivity or fixed coefficients for all the surfaces.

It would also be interesting to make the system capable of differentiating between images from different months of the year, since not all surfaces remain with the same characteristics of vegetation and soil along the different seasons. So we would be able to perform more accurate calculations, according to the date on which the original AATSR image was taken.

In this chapter, this method is applied to the European territory, but it is important to remember that it was developed always keeping in mind the need to allow it to be used worldwide in a near future. Despite there are several differences that should be taken into account, we already have the right tools to do it. So this should not be a really hard process.

Although the AATSR sensor has been used as the source of the required satellite data to apply this model, one can see that this algorithm has been developed to be applied to any other similar sensors. It would be only necessary to recalculate the parameters for the required channels.

On the other hand, as we said at the beginning of this document, this method will be used as the first step, in order to be able to use a split-window algorithm (with the purpose of correcting the atmospherical effect in satellite measures) to obtain the land surface temperature automatically.

All this effort provides an important advance in the study of climate change, weather forecasting, forest protection, fire detection, etc. All of them processes that have been traditionally tedious and cost intensive to implement, while remote sensing allow us to make this studies in a faster and more affordable way. This can be determinant in sparsely populated zones, areas of difficult access such as large forests and especially in developing countries, whose resources are often more limited.

Acknowledgments

This work was financed by the Generalitat Valenciana (project PROMETEO/2009/086, and contract of Eduardo Caselles) and the Spanish Ministerio de Ciencia e Innovación (project CGL2010-17577). AATSR data were provided by European Space Agency under Cat-1 project 3466. The authors thank ESA and the ESA GLOBCOVER Project, led by MEDIAS-France, for the GLOBCOVER classification data. The authors also wish to thank the suggestions made by Prof. C. Coll, Dr. J. M. Sánchez, Dr. M. Mira, J.M. Galve, V. García, M.M. Bisquert and C. Doña, of the University of Valencia.