1. Introduction

Up to now many Bidirectional Reflectance Distribution Functions (BRDFs) have been developed to describe the bidirectional reflectance in visible and near infrared channels as a function of both illumination and view geometries [1–6]. The semi-empirical kernel-driven models [2–5] have been proven successfully in application to AVHRR, Polarization and Directionality of Earth Reflectance (POLDER), MODIS, Multi-Angle Imaging Spectra-Radiometer (MISR), laboratory, and field-measured multi-angular reflectance data and have been shown to fit observed BRDF data well [2–6]. Only a few works focused on the BRDF modeling in the mid-infrared region (MIR), but all of them have aimed to estimate the emissivity in MIR from the bidirectional reflectance derived from AVHRR and MSG/SEVIRI data [7–9].

This paper will be devoted to estimate the land surface directional emissivity in MIR channel from the bidirectional reflectance derived from MODIS data in two adjacent MIR channels. Section 2 recalls the methodology to retrieve directional emissivity in MIR channel. Section 3 describes the study area, MODIS data and data processing. Section 4 presents some preliminary results and cross-validation with the MODIS land surface temperature/emissivity product MYD11B1 data. Finally, conclusions are given in the last section.

2. Determination of directional emissivity in MIR channel from MODIS data

2.1 Retrieval of the bidirectional reflectivity in MIR channel from MODIS data

The instrument MODIS onboard Terra and Agua satellites has two adjacent MIR channels 22 and 23 centered at 3.97 um and 4.06 um respectively. Based on the difference in the solar reflection in these two channels, and assuming that the surface bidirectional reflectivities are equal in channels 22 and 23, and that the ground brightness temperatures in these two adjacent channels are the same if the contribution of the direct solar radiation is not considered, Tang and Li [10] developed a method to retrieve the bidirectional reflectivity (ρ _b) in the MIR channel from MODIS channels 22 and 23 with

where B is the Planck function, T _{g_22} is the daytime ground brightness temperature of MODIS channel 22, R ₂₂ ^s is the solar irradiance at ground level in MODIS channel 22, T _g ⁰ is the MIR ground brightness temperature without the contribution of the solar direct beam and can be estimated from the ground brightness temperatures T _{g_22} and T _{g_23} in the channels 22 and 23 using

in which the coefficients a ₁-a ₃ are dependent only on the Solar Zenith Angle (SZA). More details concerning both the development and the application of this method with MODIS data can be found in [10].

2.2 Estimation of the directional emissivity in MIR channel from bidirectional reflectivity

For an opaque medium in thermal equilibrium, the directional emissivity ε(θ) is related to the hemispherical directional reflectance ρ _h(θ) by the Kirchhoff’s law as:

where θ is the viewing zenith angle, θ_i is the incident radiation angle, φ is the relative azimuth angle between the observation and incident directions, and ρ _b is the bidirectional reflectance of land surface in MIR channel retrieved with Eq. (1).

Based on the theory that land surface reflectance typically consists of three components: the isotropic scattering, the volumetric scattering and the geometric-optical surface scattering, a kernel-driven BRDF model, the RossThick-LiSparse-R model, was proposed to describe the non-Lambertian reflective behavior of land surface in visible and near-infrared regions [2, 4, 6]:

where k _iso is the isotropic scattering term, k _vol is the coefficient of the Roujean’s volumetric kernel f _vol, and k _geo is the coefficient of the LiSparse-R geometric kernel f _geo.

For a plane-parallel dense vegetation canopy with uniform leaf angle distribution, and equal leaf reflectance and transmittance, the Roujean’s volumetric kernel [2] is given by

where ζ is the phase angle, related to the conventional angles by

Considering the mutual shadowing between different protrusions of vegetation canopy, the reciprocal LiSparse geometric kernel f _geo derived by Wanner et al. [3] and modified by Lucht [6] is employed

where G(θ, θ_i, φ) is the overlap area between the view and solar shadows and given by

in which $\theta \text{'}={\tan }^{-1}\left(\frac{b}{r}\tan \theta \right),$, ${\theta \text{'}}_i={\tan }^{-1}\left(\frac{b}{r}\tan {\theta }_i\right)$, cosξ' = cosθ' cosθ_i' + sinθ' sinθ_i'. cosφ $t={\cos }^{-1}\left(\frac{h}{b}\frac{\sqrt{D^2+{\left(\tan \theta \text{'}\tan {\theta \text{'}}_i\sin \phi \right)}^2}}{\sec \theta \text{'}+\sec {\theta \text{'}}_i}\right)\mathrm{with}D=\sqrt{{\tan }^2\theta \text{'}+{\tan }^2{\theta \text{'}}_i-2\tan {\theta }_i^{\text{'}}\cos \phi }$ h/b and b/r are the dimensionless crown relative height and shape parameters, respectively.

According to Eqs. (3), (4) and (5) and assuming that the BRDF shapes in the MIR spectral region are the same as the ones in visible and near-infrared regions [11], the directional emissivity in MODIS MIR channel is given by

with If_x (θ) = ∫₀ ^2π f_x (θ,θ_i,φ)sin(θ_i)cos(θ_i)dθ_idφ, in which the subscription x represents vol or geo.

As shown in Eq. (10), the integrals of If_vol(θ_v) and If_geo(θ_v) over the incident radiation angle θ_i and the relative azimuth angle φ are complicated mathematical expressions and can not be analytically derived. As used in MODIS BRDF/Albedo products, taking h/b=2 and b/r=1, i.e., the spherical crowns are separated from the ground by half their diameter, Jiang and Li [11] showed numerically that the integrals of the Roujean’s volumetric kernel f _vol (Eq. 6) and the reciprocal LiSparse geometric kernel f _geo (Eq. 8) can be written with a good approximation as

It should be noted that if a series of ρ _b with different angular configurations are retrieved from the MODIS data using Eq. (1), one can get the parameters k _iso, k _vol and k _geo from Eq. (5). Knowing these three parameters, the directional emissivity in MIR channel can be obtained with Eqs. (10), (11) and (12).

3. Study area and data processing

A region of Egypt and Israel with latitude from 28.0° N to 32.0° N and longitude from 30.0° E to 36.0° E was chosen in this study. Figure 1 shows the land use map of this study area generated from MODIS land cover type 2004 L3 global 1 km product MOD12Q1 (http://edcdaac.usgs.gov/modis/mod12q1v4.asp) and classified by the International Geosphere-Biosphere Programme (IGBP). From this figure, we can see that the major land cover types in this area are barren or sparsely vegetated, croplands, and open shrubland. Since the classification scheme of IGBP does not include bare soil surface, to discriminate from barren or sparsely vegetated, this surface type will be used in our following work. The reason to choose this region is that a series of cloud-free MODIS data are available over this region from July 12 to July 30 of 2005. Therefore, the three parameters k _iso, k _vol and k _geo in Eq. (5) can be determined with the retrieved bidirectional reflectivity ρ _b under the assumption that the land surface remains unchanged during this period.

The MYD021KM, MYD03 and MYD35_L2 product files provided by the NASA Goddard Space Flight Center (GSFC) Level 1 and Atmosphere Archive and Distribution System (LAADS) (http://ladsweb.nascom.nasa.gov/data/) were used in our work. The MYD021KM data, calibrated Earth View data at 1 km resolution by the MODIS Characterization and Support Team (MCST), are the Top of Atmosphere (TOA) radiances and reflectances. The geolocation dataset, MYD03, provides latitude, longitude, ground elevation, solar zenith and azimuth angles, and satellite zenith and azimuth angles for each 1 km sample. The MYD35_L2 is a cloud mask product which gives a clear-sky confidence level (clear, probably clear, uncertain, cloudy) to each IFOV. More details about these product files can be found in [12]. Ten days MODIS data with cloud-free conditions at the moment of MODIS overpasses, from July 12 to July 30 of 2005 were selected. Table 1 gives the dates and acquisition times of these ten days MODIS data. The European Centre of Median-range Weather Forecast (ECMWF) reanalysis (ERA) operational deterministic model data directly obtained from the French Meteorological Center with latticed resolution of 0.5° in both latitude and longitude [13] were used to perform atmospheric corrections for MODIS MIR data in this work. In addition, taking into account the real atmospheric path length between the surface and the satellite, global DEM data at 30 arc-s (1 km) resolutions (http://edc.usgs.gov/products/elevation/gtopo30/gtopo30.html) were also used.

 

Fig. 1. Land use map of the study area generated from MODIS land cover type 2004 L3 global 1 km product (MOD12Q1) and classified by IGBP classification scheme.

Download Full Size | PDF

Table 1. Date and acquisition time for ten days MODIS data used in this study

View Table | View all tables in this article

Based on the clear-sky confidence level (clear, probably clear, uncertain, cloudy) assigned to each IFOV in MYD35_L2, clear and probably clear pixels were taken as clear, and uncertain and cloudy pixels were taken as cloudy in our study. The cloudy pixels assigned in this study were then firstly screened out in the retrieval of ρ _b. Since the satellite instrument measures only the radiances at the Top of Atmosphere (TOA), the data acquired by MODIS MIR channels 22 and 23 have to be corrected for the atmospheric effects in order to obtain the radiances or brightness temperatures at ground level. These atmospheric corrections were performed using the atmospheric radiative transfer model-MODTRAN 4 with the ECMWF data and DEM data. Selection of ECMWF output data as atmospheric profiles is due to the fact that the MIR channels 22 and 23 are not too sensitive to the change of water vapor content in the atmosphere. More details of atmospheric corrections for MODIS MIR channels can be found in [10]. After having performed the atmospheric corrections, the bidirectional reflectances in MODIS MIR channel 22 can be estimated with Eqs. (1) and (2).

The three parameters k _iso, k _vol and k _geo in Eq. (5) for each pixel are then determined by a Levenberg-Marquardt minimization scheme with the retrieved ρ _b and corresponding illumination and view angles extracted from MYD03 data. Finally, the directional emissivities at each view zenith angle for MODIS MIR channel are obtained with Eq. (10).

4. Results and validations

The objective of the present work is to estimate the directional emissivity from MODIS MIR channels. Figure 2 gives an example of the retrieved directional emissivity map for July 24, 2005 at 11:25 UTC time.

 

Fig. 2. Map of the directional emissivity in MIR channel for July 24, 2005.

Download Full Size | PDF

As shown in Fig. 1, points A, B, C, and D marked in Fig. 2 represent bare soil, open shrubland, barren or sparsely vegetated, and croplands surfaces respectively. For the entire study area, the directional emissivities in MODIS MIR channel 22 vary from 0.6 to 1.0, and they are usually less than 0.80 over the bare areas, while the opposite is observed over the vegetated areas. Figure 3 illustrates histograms of the estimated directional emissivities for the four major land covers (bare soil, open shrubland, barren or sparsely vegetated, and croplands surfaces) in the entire study area. As displayed in Fig. 3, the directional emissivity in MIR channel varies from 0.67 to 0.78 with mean=0.73 and standard deviation (std)=0.021 for the bare soil surfaces, and from 0.83 to 0.94 with mean=0.89 and std=0.019 for the open shrubland, while for barren or sparsely vegetated surfaces the directional emissivity in MIR channel ranges from 0.89 to 0.98 with mean value of 0.94 and std=0.013, and for croplands, the emissivity in MIR channel is the highest and ranges from 0.92 to 0.99 with mean=0.97 and std=0.012.

 

Fig. 3. Histogram of the directional emissivity in MIR channel estimated from MODIS data for the major land cover types in the study area. Std=standard deviation.

Download Full Size | PDF

Figure 4 displays the sun and satellite zenith and azimuth angles in polar representation at four locations for ten clear days from July 12 to July 30 of 2005 which we used to retrieve the directional emissivity. For these locations, the sun is in the West direction and coincides nearly with the satellite along track direction. The observation directions are almost in the principal plane and lie in the east and west directions according to the instrument scanning directions. It should be pointed out here that although the change of solar zenith angle is very small during this period for a given pixel, the viewing zenith angle of each pixel (location) varies significantly from 0° to 60° from July 12 to July 30. We can, consequently, get the parameters k _iso, k _vol and k _geo using Eq. (5) with a series of ρ _b and different angular configurations.

Figure 5(a) shows the comparison of the bidirectional reflectance ρ _b estimated directly from MODIS MIR data (Eq. 1) with ρ _b modeled using Eq. (5) at four locations for the ten clear days. The Root Mean Square Error (RMSE) and Mean Error (ME) are respectively 0.005 and zero. From this figure one can notice that the bidirectional reflectances for locations A and B are larger, while for locations C and D, they are relatively smaller. In addition, Fig. 5(b) gives the histogram of the differences between the retrieved and modeled bidirectional reflectances for the entire study area. From this figure, one can see that the difference of the retrieved and modeled ρ _b ranges from -0.03 to 0.03 with mean of -0.001 and standard deviation of 0.008.

 

Fig. 4. Sun and satellite zenith and azimuth angles in polar representation at four locations for ten clear days during the period of July 12 to July 30 of 2005.

Download Full Size | PDF

 

Fig. 5. Comparison of the bidirectional reflectances estimated using Eq. (1) with those modeled using Eq. (5): (a) for locations A, B, C and D, (b) for the entire area.

Download Full Size | PDF

Table 2 gives the fitting parameters of k _iso, k _vol, and k _geo in Eq. (5) for locations A, B, C and D. In addition, values of the Normalized Difference Vegetation Index (NDVI) derived with the TOA reflectances in near-infrared and red channels from MYD021KM data for these four locations on July 24, 2005, are also given in this table.

Table 2. Fitting parameters k _iso, k _vol, and k _geo in Eq. (5) for locations A, B, C, and D

View Table | View all tables in this article

To preliminarily validate the directional emissivity estimated using the present method, the MODIS land surface temperature/emissivity product MYD11B1 data were used in our investigation. Taking into account that the MYD11B1 product provides the land surface emissivity values at 5 km resolution, the mean value of estimated directional emissivities for 5×5 pixels with 1 km resolution was selected to match the one from MYD11B1 data with regard to the nearest latitude and longitude coordinates. Figure 6(a) displays the directional emissivities estimated using the present method versus those extracted from MODIS land surface temperature/emissivity product MYD11B1 data at four locations for the ten clear days. From this figure we can see that the Mean Error (ME) and the Root Mean Square Error (RMSE) are 0.002 and 0.021 respectively. In addition, Fig. 6(b) shows cross-comparisons of the estimated directional emissivities and those extracted from MYD11B1 data for the entire region for these ten clear days. The ME and RMSE between the directional emissivities estimated in this study and those extracted from MYD11B1 data are of -0.007 and 0.024 respectively. The result of this comparison shows that, at least for our cases, the method described in this paper for estimating the directional emissivity in MIR channel gives results comparable to those in MYD11B1 product.

 

Fig. 6. Comparison of the directional emissivities estimated from MODIS MIR channels using Eqs. (1), (5) and (10) with those from MYD11B1 product for ten clear days during July 12 to 30, 2005: (a) for four locations (b) for entire regions.

Download Full Size | PDF

5. Conclusions

In this work, the directional emissivity in MODIS MIR channel has been estimated with the retrieved bidirectional reflectivity and the RossThick-LiSparse-R model. Ten days of MODIS MIR data with cloud-free conditions at the moment of MODIS overpasses, from July 12 to July 30 of 2005 were used to determine the three parameters k _iso, k _vol and k _geo in Eq. (5) for each pixel. The directional emissivities of these days were mapped for a region of Egypt and Israel with latitude varying from 28.0° N to 32.0° N and longitude from 30.0° E to 36.0° E.

In order to show the retrieval accuracy of the proposed method, the MODIS land surface temperature/emissivity product MYD11B1 data have been used to cross-validate preliminarily the directional emissivities derived directly from MODIS MIR data with the method presented in this paper. The results of this comparison showed that, at least for our cases, the proposed method for estimating the directional emissivity gives results comparable to those of MYD11B1 product with Mean Error =-0.007 and Root Mean Square Error =0.024.