Demonstration of uncertainty resulting from MODerate resolution Imaging Spectroradiometer-like geometries: an albedo case study

ABSTRACT

This research examines uncertainty in MODerate resolution Imaging Spectroradiometer (MODIS) observations, and demonstrates the direct influence of geometric distortions resulting from the standard practice of geolocating swath observations. MODIS observations vary dependent on the ground sample distance, which varies dependent on the view zenith angle that changes with every orbit. MODIS Level 2G (L2G) land products are generated by applying a geolocation algorithm that resamples the variable observation geometries to a consistent grid of fixed pixel size and location, a process which itself introduces variability associated with the changing observational footprint. For this study, broadband albedo was simulated for five validation sites, representing five distinct land cover types, exhibiting quantifiable variability, with additional seasonal variability exhibited in some sites. All site simulations exhibit compounded uncertainty attributable to the geometric distortion sufficient to influence climate models (i.e. ranges from 0.01 to 0.045 albedo). These results indicate there is a minimum level of uncertainty associated with the variable geometry that should be factored into L2G-based products, particularly for nominal 250 m band data. Aggregating the data to coarser resolutions and smoothing the data through average resampling can mitigate the uncertainty.

KEYWORDS:

• Earth observation
• geography
• land cover
• remote sensing
• visible/infrared remote sensing

1. Introduction

Coarse resolution sensors (>250 m) commonly used for land remote sensing, such as the NOAA AVHRR, SPOT VGT, and MODerate resolution Imaging Spectroradiometer (MODIS), have broad swath-widths (2200–2600 km). The resulting off-nadir viewing allows for near-daily viewing, increasing the opportunity for cloud-free observations (Townshend and Justice 2002). However, the off-nadir viewing results in variable view geometries, which need to be addressed in data pre-processing and production, and provides an inherent uncertainty in the resulting products. This paper examines the effects, and potential significance to climate models, of measurement uncertainty imposed by the variable viewing geometry of the MODIS instrument using a hypothetical albedo product. For the purposes of this study, albedo is used as a convenient and familiar physical quantity which is well established within the comprehension and application of the scientific community.

Albedo is the fraction of incident solar irradiance reflected by the Earth's surface across the solar spectrum (Goward 2005), measured as a ratio of upwelling to downwelling radiance, providing estimates of the energy budget available for driving heat fluxes in climate system modeling. Radiative forcing is the change imposed to Earth's radiation balance as a result of some perturbation of the surface (Hansen, Sato, and Ruedy 1997), measured in Watts per meter squared (Wm⁻²), providing a useful measure for positive (warming) and negative (cooling) changes, and can be derived from measurements of albedo. The relationship of albedo in a basic net radiation budget model is:(1) where is the surface albedo, is directly absorbed solar radiation, is surface emissivity, is downwelling long-wave atmospheric emission ( = air temperature), and is upwelling long-wave surface emission ( = surface temperature).

Monteith and Unsworth (2008) demonstrate the linear relationship between irradiance and air temperature for both direct and diffuse radiative scattering models; multiple studies have either included or modeled their climate system response to changing land albedos with respect to direct temperature effects and carbon sequestration (Dixon et al. 2003; Montenegro et al. 2009; Schurgers et al. 2008; Waggoner and Sokolik 2010).

The role and importance of albedo in Earth Systems Science have been established since the early days of remote sensing (Sagan, Toon, and Pollack 1979), and are now key to the ongoing discussion regarding climate change and anthropogenic effects upon Earth systems (Goward 2005; Lacis et al. 2010). Recorded albedo estimates derived from satellite observations contribute to our understanding of changes in the climate, through a relation to the net radiative flux over different landscapes. The absolute error for satisfactory albedo parameters in climate models ranges from to (Oleson et al. 2003). As climate models continue to improve their spatial and temporal resolution and incorporate new data sources (IPCC 2007), identifying and minimizing model error become increasingly necessary.

Two studies, each using simulated cropland land cover class albedos, exemplify the dependency between albedo accuracy and model uncertainty. Myhre and Myhre (2003) ran sequential models using albedo values from 0.15, 0.18, and 0.20 to simulate baseline albedo scenarios representative of the potential natural vegetation (PNV); the resultant radiative forcing was calculated at , and Wm⁻² cooling, respectively. Matthews et al. (2003) ran models using albedo values between 0.17 and 0.20 to simulate historical baseline land covers from the 1700s CE to 1992, demonstrating increased radiative forcing from to Wm⁻², calculated as a difference in cooling of −0.09°C and −0.17°C, respectively. Both studies conclude that model outcomes are dependent on small changes in albedo parameters and describe a general relationship where albedo corresponds to and −0.03°C cooling.

The MODIS, flown on two polar orbiting platforms, has provided useful information regarding planetary albedo since the early 2000s (Schaaf et al. 2002). The Fifth Edition of the Intergovernmental Panel on Climate Change (IPCC) describes the uncertainties of the albedo parameters used in global climate models (Myhre et al. 2013), exhibiting the sensitivity of model runs to changes in albedo parameters, and demonstrating the dependency of model outcomes on the accuracy of albedo parameters, some of which are derived from MODIS. Calculations of albedo are sensitive to observational characteristics, such as atmosphere and land cover. For example, the Myhre (2009) albedo study demonstrates diminished negative forcing, to Wm⁻², when atmospheric models incorporate aerosol optical depth (AOD) from MODIS observations. Liang, Strahler, and Walthall (1999) report albedo uncertainty exceeding attributed to sources including changing AOD conditions, variable zenith angles, and complex landscapes.

Given the sensitivity of climate models to small changes in albedo, this paper simulates a series of MODIS-like gridded broadband albedo for several land cover types to test the measurement consistency of a given gridded pixel through changing view geometries. In the theoretical space-time designed for this study, only the MODIS satellite view angle changes relative to the Earth's surface, while the observed relationship between the Earth's surface (study sites) and all other potential effects (angles ofincidence, temperature, atmosphere, etc.) to the calculation of albedo are assumed constant. This paper demonstrates the level of measurement uncertainty imposed by the variable view geometry of the MODIS instrument and geolocation process at different spatial scales using broadband albedo as a measure. It is important to note that this paper is not in any regard an examination or validation of existing albedo products or their derivation protocols. Therefore, the simulation of the hypothetical albedo does not employ the same production methods or quality assessment as those products. Reference to any given published product accuracies is solely to compare the relative error contribution.

2. MODIS observation characteristics

MODIS is a whisk-broom system with a scanning mirror providing a wide observation swath (2330 km) for 36 spectral bands at nominal detector pixel resolutions of 250, 500, and 1000 m (Barnes, Pagano, and Salomonson 1998). The wide field-of-view (±66°) and high-temporal frequency (at least 2 passes a day at the equator) of MODIS provides near-daily global coverage (Salomonson et al. 2000; Xiong and Barnes 2006), but comes at the cost of a loss in signal accuracy due to known effects of pixel size growth at increasing view zenith angles (VZAs). VZA is the angle formed between a line perpendicular to the observed surface (local zenith) and a line pointing back to the instrument.

The ground-projected sample interval (the spacing between pixels), known as the ground sample distance (GSD), is a common measure of instrument spatial resolution (Schowengerdt 2007). The reported values of 250, 500, and 1000 m correspond to the nominal GSD of the MODIS instrument at nadir and do not reflect the continuously varying VZA, the overall effect of which is a growth in the GSD of approximately 2.0 times along-track and 4.8 times along-scan at the maximum observation angle. The pixel growth ultimately leads to overlapping scans toward the edge of scan, known as the bow-tie effect (Nishihama et al. 1997). The constant scanning of the mirror forms a triangular line spread function (LSF) (a sensor's one-dimensional spatial response to the measured physical signal) integrating radiometric signal from the immediate preceding and subsequent neighbor pixels (12.5%) to the center pixel (75%) along-scan (Nishihama et al. 1997; Wolfe et al. 2002).

To help facilitate use, all higher order MODIS products are assigned to a geo-referenced grid using a custom sinusoidal projection with a fixed pixel size for each of the available nominal resolutions. MODIS geolocation assumes a fixed spatial size based on the subdivision of a 30-arc second grid providing a seamless global grid with resolutions of , , and (Wolfe et al. 2002). Observations made at resolution (e.g. 250 m) are assigned to a co-located () coordinate in the geolocation grid of the same resolution. Observational coverage (obscov) is the integrated signal response during a sample interval overlapping a grid cell divided by the full observation footprint, which accounts for the changing sample sizes at increasing GSD. Grid pixels are assigned the value of those observations (swath pixels) with the highest obscov values (Wolfe, Roy, and Vermote 1998; Yang and Wolfe 2001). However, the reported mean obscov for MODIS grid data is less than 30% (Tan et al. 2006), meaning the value assigned to any given grid pixel is only a portion of the spatial signal from the larger observation pixel, and not necessarily representative of the small geographic area defined by the grid location. Further, single observations are often assigned to more than one grid pixel at large VZAs.

There is no specific problem with the assignment of single or partial values to multiple grid cells. Conversely, there is no benefit from producing a finer resolution product, given the cell will only ever represent a portion of spatial signal present because the GSD of the instrument is not consistent and is always greater than the grid cell size. The portion of the spatial signal will also differ as a function of location, since the observation pixel center location and associated geometry will be different for each orbit between nadir-to-nadir views, making the defined grid space an arbitrary subdivision of the swath data.

The results from Campagnolo and Montaño (2014) introduced an inconsistent, location-dependent effective GSD that is far different from the nominal grid resolutions suggested by product descriptions. The author concludes that the effective resolution of MODIS 250 m grid data varies between 344 and 835 m along-scan and behave in the same basic manner as swath data (i.e. the effective resolution grows in the same fashion as swath data).

Many MODIS products are generated as discrete temporal composites, using only the single best quality per-pixel observation (e.g. cloud-free and/or minimum VZA) for a given composite period – usually 8–16 days. Data in these composites ultimately suffer from geometric distortion in that a seemingly spatially continuous scene nonetheless represents a spatially discontinuous surface comprising observations from any number of disparate GSDs.

Figure 1 illustrates the different days of data and the associated VZAs used in a typical discrete composite MOD13 NDVI product. The composite day of year seen in the figure is an interesting illustration of the false assumption of spatial continuity of data within a composite scene. There are 8 days worth of data in the scene that are not necessarily spatially contiguous. Further, and of particular import to this study, the range of VZAs in the composite is far larger than one might expect. Specifically, the MOD13 product series applies a minimum VZA criterion in the selection of best quality pixels (Solano et al. 2010). However, there are clearly circumstances overriding these criteria (e.g. a limited number of cloud-free observations) making the use of less desirable observation days (e.g. those with large VZAs) a necessity.

Figure 1. Sample illustration of days used in compositing the MODIS MYD13Q1 NDVI product from the 27 April 2008 (day 116) 16-day composite. The top figure shows the VZAs associated with the days used. The lower figure shows the day within the period that was used in the composite. The MOD13 product uses a minimum VZA as part of the quality selection, but nonetheless some large angles are used as a result of coincident factors such as persistent cloud contamination.

3. Data and methods

Assessing the impact of geometric artifacts on measures of albedo requires isolation of the geometric influence from other potential sources of error. To accomplish this, synthetic time series of MODIS data were generated from 11 individual scenes of Landsat 5 TM images for five earth observation system (EOS) validation sites (Morisette, Privette, and Justice 2002), representing five distinct land cover types (Table 1).

Table 1. EOS validation sites and landsat 5 TM data used in the study.

Site Name	Lat	Lon	Path/row	Date	Land cover
Bondville, IL	40.00	−88.29	22/32	2010-04-14	Broadleaf Cropland
				2010-07-03
				2010-11-08
Harvard Forest LTER, MA	42.54	−72.18	13/30	2008-04-26	Broadleaf Forest
				2008-06-12
				2006-10-13
Konza Prairie LTER, KS	39.08	−96.56	28/33	2011-03-10	Grassland/Cereal Crop
				2011-07-01
				2011-09-02
Jornada LTER, NM	32.60	−106.86	33/37	2011-04-30	Shrubland/Woodland
Lake Tahoe, NV	39.17	−120.104	43/33	2009-09-21	Needleleaf Forest/Clear Water

Hypothetical broadband albedo was calculated using the simple linear model and methods for non-snow land cover described in Brest and Goward (1987). The hypothetical broadband albedo method uses a linear combination of visible and near-infrared (NIR) bands calibrated to top-of-atmosphere reflectance to capture the characteristic absorption of visible and reflectance of NIR spectrum demonstrated by vegetation. Brest and Goward (1987) considered either Landsat MSS bands 4 or 5 for the visible portion of the spectrum, but ultimately chose band 4 as the most representative reflectance in the visible spectrum between 0.3 and 0.725 μm. The more precise spectral response in the visible/NIR border region in the Landsat 5 TM sensor, along with the coincident bandwidths in 250 m MODIS bands 1 and 2, led to the author's decision to substitute the original Landsat MSS band 4 with Landsat TM band 3.

The simple broadband model is sufficient to serve as a context from which geometric variability can be assessed. To the degree possible, all other variables (e.g. changing angle of incidence, and atmosphere) have been constrained; variables not contributory to the generation of synthetic albedo for the purposes of this study have been assumed as constant.

The method first uses a simple ratio (SR) of band 4 divided by band 3, where a ratio greater than 2 indicates vegetation versus non-vegetation. The broadband albedo is then calculated using the appropriate equation described in Equation (2).(2)

The surface composition and spatial variability of the land cover at any site will influence the amount of spatial signal incorporated into any given GSD; a uniform, homogeneous canopy will present a more stable spatial signal over coarsening scales than a highly variable, heterogeneous landscape with small crop fields. Spatial variation further depends on seasonality for those canopies that experience different levels of canopy green-up and senescence. Therefore, seasonal profiles were generated for those sites expected to demonstrate seasonally variable phenology: Harvard Forest (broadleaf forest), Konza Prairie (grassland/cereal crop), and Bondville (broadleaf cropland). Images containing snow cover are beyond the scope of this analysis, and were excluded.

3.1. Study areas

The Bondville, IL, EOS agricultural site is located in the Midwest USA and is characterized by agricultural land with generally rectangular fields of agricultural land that can be in various states of development based on both time of year and crop type in each field. A high level of spatial heterogeneity can be expected from this site as a function of the GSD being larger than the field size from day to day.

The Landsat albedo pixels centered around the EOS validation site extending in a circle to 250, 500, 1000, and 2000 m (Figure 2(a)) and for each season were extracted to understand the level of spatial heterogeneity of the study sites at varying scales. The resultant distributions for each scale in Figure 2 are used to understand the level of spatial heterogeneity of the study sites. Non-normal distributions will incorporate varying amounts of outlier data as the MODIS GSD changes, causing irregularity in the simulated MODIS albedo; a normal distribution (i.e. spatially homogeneous) provides a more stable input signal for each simulated footprint, thus less variability as the GSD changes. The distributions for all seasons at the Bondville site indicate an irregular distribution at finer scales that tend to smooth at coarser resolutions, indicating inherent spatial heterogeneity at the site.

Figure 2. Distributions of the source Landsat albedo for each season at increasing spatial scales (250–2000 m) centered around the site.

The Harvard Forest LTER (EOS Test Site) is a deciduous forest in the northeast USA, which would be expected to demonstrate some seasonality as a result of leaf-on/leaf-off periods. The histograms for the three seasons in Figure 3 are generally more normally distributed than for the Bondville site, suggesting this site is more spatially homogeneous. However, the distributions demonstrate a skew at finer resolutions (250 and 500 m), particularly in the Summer period, which can manifest as greater variability with varying GSD.

Figure 3. Histograms of albedo frequency for the different seasons and scales for the Harvard Forest site.

The Konza Prairie LTER site is a grassland land cover class in the Midwest USA. The histogram distributions in Figure 4 have a relatively normal distribution for all scales, similar to that of Harvard Forest. The Summer period shows a skew at finer resolution, though the frequency is denser than the skew in the Harvard site, which suggests the site will continue to be homogeneous.

Figure 4. Histograms of albedo frequency for the different seasons and scales for the Konza Prairie site.

Land cover at the Lake Tahoe site was not expected to demonstrate high levels of seasonality. Therefore, only a single Fall image was simulated for Lake Tahoe. The Lake Tahoe EOS site is a clear water and needle leaf land cover class in the western USA. The Lake Tahoe site was chosen for the interesting feature of having the exact site location defined in the lake, but the aggregation of pixels ultimately incorporating land features at coarser resolutions as illustrated in Figure 5(a). The Jornada LTER, located in the southwest USA, is classified as a semi-arid grass/scrubland. A single Summer image was selected for the Jornada site, since seasonality was not expected to influence the land cover.

Figure 5. (a) for Lake Tahoe illustrates the composition of pixels at coarser resolutions and the level of land features that are added to the signal. The relative homogeneity of the Jornada site can be seen in (c).

The histograms for both Lake Tahoe and Jornada are seen in Figures 5(b) and 5(d), respectively. The histograms for the Lake Tahoe site are very irregular, which is an artifact of the incorporation of very different land albedos within the scene at coarser scales, and is highly spatially homogeneous at fine scales. The spatial homogeneity of the Jornada site is indicated through the generally normal distributions of albedo at all scales.

3.2. MODIS albedo simulation

MODIS bands 1 and 2 used for the simulation are nominal 250 m bands; thus, the decision was made to start the analysis for this research at the finest available pixel size. The process of generating the synthetic MODIS data from the Landsat scenes incorporates the known spatial integration of signal by creating a triangular weighting scheme for all Landsat pixels that fall into a given MODIS swath observation GSD. The simulation is accomplished in several steps.

A MODIS observation center for the 250 m resolution must first be determined from the input geolocation file (MOD03/MYD03), which is stored at 1 km resolution. The cell centers for the first along-track 1 km earth scans are co-registered with the cell centers of the second 500 m and fourth 250 m spatial samples. Figure 6 from Nishihama et al. (1997, Figures 2–11) shows the relationship between the spatial resolutions. Accurate geolocation must occur for each scan independently due to the overlap of successive scans. Therefore, each scan is calculated individually by dividing the number of rows in the MOD03/MYD03 file by 10, that is, the number of scans in each file. The 250 m longitude and latitude for each pixel center were then calculated using bilinear interpolation of the 1 km longitude and latitude values. The first two and last two rows in the 250 m geolocation fall outside the actual 1 km frame, so these values were approximated using bilinear extrapolation.

Figure 6. The spatial relationship between 1 km geolocation data and the 500 and 250 m resolutions. The cell center of the 1 km data is co-registered with the second and fourth 500 and 250 m pixels, respectively. An accurate geolocation for each of the resolutions can be determined though bilinear interpolation of the 1 km values. (This figure is derived from Nishihama et al. (1997, Figures 2–11).)

For this study, the observation footprint of each discrete pixel was derived from the 250 m cell centers per the methods described in Tan et al. (2006, Appendix B). The cell centers themselves were derived from the specific viewing geometries of the MOD03/MYD03 geolocation file (i.e. approximate scan angle, sensor azimuth, and known sensor altitude). The variability of the viewing geometries result in variability amongst the discrete pixel sizes.

Radiometric signal value was assigned to each discrete pixel, in a simulation of the actual functioning of the triangular along-scan LSF of the MODIS instrument (Nishihama et al. 1997). Both the simulated signal value and the actual signal value are aggregates, not discrete values, composed of signals from the pixel of interest (75%) and both neighboring pixels along-scan (12.5%). In this study, a weighting scheme of W = [0.2, 0.4, 0.6, 0.8, 1.0, 0.8, 0.6, 0.4, 0.2] was applied to approximate the relative contribution of signal to each discrete pixel along-scan, while equal weights were assigned to each discrete pixel along-track, per established regularity (Barnes, Pagano, and Salomonson 1998), as illustrated in Figure 7.

Figure 7. Illustration of swath observation simulation. Landsat pixels are assigned weights from a 5 × 9 grid centered inside the swath pixel and incorporating signal from neighboring pixels. The value for a given pixel is then the weighted average of all Landsat pixels. The simulated swath observation is consistent with the LSF of the MODIS instrument.

This matrix of weighted values along-scan (weighted columns) and regular values along-track (regular rows) resulted in a semi-weighted 5-row × 9-column grid. In order to represent the Landsat pixels within the simulated MODIS context, Landsat pixels found to be located inside the simulated MODIS observation footprint of each discrete pixel, and adjacent simulated MODIS pixels, were assigned weights from the nearest weighted column; nearest neighbor resampling allowed for computational simplicity. Along-track weights being accepted as regular were treated as constant and discrete; for computational simplicity no effort was made to interpolate between column weights. The final observation value was calculated using the weighted average of all Landsat pixels within the semi-weighted grid, in approximation of the established 12.5%/75%/12.5% neighborhood contribution and the triangular shape of the MODIS LSF (Nishihama et al. 1997).

The individual Landsat albedo scenes were transformed to MODIS swath geometries for each MOD03/MYD03 file intersecting a given site during an -day period, representing at least two nadir-to-nadir periods (see Table 2). The resultant synthetic MODIS swath data were maintained in a WGS84 geographic projection, preserving the variable pixel size for the daily swath. The synthetic swath data were projected from geographic to the custom MODIS sinusoidal projection with a pixel resolution of 250 m. The choice was made to keep the grid resolution at an even 250 m for easy manipulation even though the actual MODIS grid resolution is . The reprojection uses a nearest neighbor pixel assignment method to approximate the obscov method used in actual MODIS product generation. The assignment difference between the resampling methods is expected to be small and is consistent with tools such as the LPDAAC MODIS Reprojection Tool Swath (LPDAAC User Services 2015).

Table 2. MOD/MYD03 time periods used for simulating data.

Site	Platform	Season	Begin Date	End Date
Bondville	Aqua	Fall	2010-11-08	2010-12-09
	Aqua	Summer	2010-07-03	2010-08-03
	Aqua	Spring	2010-04-14	2010-05-14
Harvard Forest	Aqua	Fall	2006-10-13	2006-11-14
	Aqua	Summer	2010-07-03	2010-08-03
	Aqua	Spring	2010-04-14	2010-05-14
Konza Prairie	Terra	Fall	2011-09-02	2011-10-01
	Terra	Summer	2011-06-30	2011-07-31
	Terra	Spring	2011-03-10	2011-04-08
Jornada LTER	Terra	Spring	2011-04-30	2011-07-08
Lake Tahoe	Aqua	Fall	2009-09-21	2009-10-26

The daily MODIS 250 m Surface Reflectance (MOD09GQ), used as input for many higher order products including albedo, processing steam selects the highest quality individual pixels from multiple orbits for the first layer based on quality criteria including overall observation quality, maximum obscov, and minimum VZA (Vermote, Kotchenova, and Ray 2011). However, the simulated data in this study were considered to have equal observation quality and obscov was not generated. Therefore, only the minimum VZA criteria was applied to select individual daily observations from multiple orbits.

Nominal 250 m band MODIS data used in coarser resolution gridded science products (e.g. MODIS Albedo MCD43 at 500 m resolution) are aggregated using the weighted average of all observations intersecting a given coarse resolution cell determined by the obscov values for the input pixels. For this study, since obscov was not generated, a simpler method was applied assuming any given resolution is an equal subdivision of the original sinusoidal grid – a 250 m pixel in the sinusoidal grid is the 4 × 4 subdivision of the 1 km grid, for example. Each resolution was calculated from the original 250 m grid data, with the simple mean value calculated from all pixels falling into the new, coarser pixel. The full process is illustrated in Figure 8.

Figure 8. The full process of synthesizing data. Clockwise from top left: (a) The original Landsat synthetic albedo at 30 m resolution. (b) The next step transforms the Landsat at 30 m resolution to daily MODIS observation geometry derived from daily MOD03 geolocation data and projected to a geographic coordinate system preserving the variable swath pixel sizes. This example uses data from Aqua and has a VZA of 26.47°, resulting in an pixel. (c) Finally, the daily synthetic swath data are gridded to the MODIS sinusoidal grid at 250 m resolution.

Although the process is not completely consistent with the production code used for MODIS land products, the process does provide an approximation of the variability. The results of the simple method demonstrate a similar effect of data smoothing through the averaging of data to coarser resolution even though the method of aggregation differs.

The final step in the process is to extract the pixel value for the pixel containing the EOS site. The sites themselves are point locations, but the albedo value was extracted from the simulated MODIS pixel containing the site. Albedo values were extracted from each pixel whose center is closest to the site coordinates using a nearest neighbor approach for each gridded scene in the series.

Two measures were used to determine variability in the data: range and standard deviation. Standard deviation provides the basic metric of variability around the observed mean values. The range of measured albedo () was a useful measure in this case to assess the geometric error. The range of signal variation from peak to trough through time is related to the changing observation GSD size as a function of the variable view geometry (VZA), since the source Landsat albedo only varies due to relative contribution in the changing GSD. Observations made at increasing VZAs will necessarily incorporate more signal from adjacent space as the observation footprint grows to be greater than the grid pixel size. The range of values for homogeneous land covers should be small, since the underlying albedo would vary little with changing GSD. Conversely, spatially heterogeneous land covers were expected to demonstrate increased ranges, since the underlying albedo will vary in response to the more complex landscape as the GSD varies in size.

4. Results

The goal of the simulation process is to hold albedo constant, with the geometric effects from the MODIS instrument as the only variable; geometric variation and inconsistent swath-to-grid assignment from orbit-to-orbit lead to differing measurements of the same underlying albedo at the same location, as the GSD changes. Data from the MODIS instruments on both Aqua and Terra are used to illustrate any differences to be seen in the viewing geometries from a morning descending orbit (Terra) versus an afternoon ascending orbit (Aqua).

An example of albedo and VZA pixel values extracted from simulated swath data at the Harvard Forest, MA site in Figure 9 clearly shows that the extreme values are not observed at the extreme range of VZAs. Rather, the extreme values are observed in the middle range of VZAs between ≈ 5° and 35°, while the smallest is near nadir. The range incorporating these values is relatively large at 0.038, though the standard deviation is far smaller at 0.009. Excluding the extreme values still retains a range of values . These measurements are made with data prior to the grid process, suggesting a significant variability from the GSD alone.

Figure 9. Simulated swath pixel albedo plotted as a function of VZA for Harvard Forest site. The results show the largest values are observed in the mid-range of VZAs between ≈ 15° and 35°. Overall variability for these simulated swath data appears low, reporting a measured standard deviation of 0.009. The range of values is larger at 0.038, which is a significant contribution to uncertainty from the variable GSD before gridding occurs.

The results for all sites and seasons during the entire simulated period (Figure 10) are organized in ascending order of increasing pixel size ranging from 250 to 2000 m in 500 m increments, showing the range and standard deviation of daily measured albedo. The results clearly demonstrate a measurable level of variability from the geometric properties of the MODIS instrument, which would carry forward to the final grid product. Reduction in both standard deviation and range is achieved through aggregating data to coarser grid resolutions, but full removal is not achieved at any resolution in the simulation. The increasing values for the Lake Tahoe site are an artifact of the center of the EOS site being located offshore, but ultimately incorporating land features at increasing grid resolutions.

Figure 10. Albedo ranges and standard deviation for all sites. The plot is arranged in ascending order of increasing pixel size from 250 to 2000 m, illustrating the albedo variability for the three seasons tested.

The overall low values for standard deviation, generally albedo, suggest that a given temporal series of data is relatively stable. However, the relationship between climate model and albedo parameters was previously described as albedo, corresponding to and −0.03°C cooling. The level of variability seen at the finest resolution, 250 m, are sufficient to influence radiative forcing in climate models attributed to the variable observation geometry of the instrument and geolocation, independent of atmospheric or albedo model-related sources.

The reduction in both range and standard deviation from aggregating the data to 500 m in this study suggests that products such as the 500 m MODIS albedo product (Schaaf et al. 2002, 2014), MCD43, incorporating multiple days of observations from both Aqua and Terra, and therefore GSD, to approximate a nadir adjusted value are not impacted greatly from the varying geometry. However, the higher values for the range indicate that significant levels of uncertainty remain when individual geometries are integrated into temporally discrete composite data products such as the NDVI in Figure 1. This applies to both assignment of varying GSD in a given scene, or from composite period to period. The variability observed from fine-resolution GSD in a scene could affect land cover classification accuracy, for example.

The focus of the above research was to demonstrate the impact of geometric distortions as a result of off-nadir viewing. The analysis shows the impact, and that, in general, the smoothing effect of aggregating the data to coarser resolutions is expected to reduce the resultant measurement variability. The analysis also revealed a relationship between land cover spatial heterogeneity, the variability in the signal, and the level of spatial aggregation required to reduce albedo variability at different sites, which is described individually in the following sections.

4.1. Bondville, IL

The results for the Bondville site (Figure 11) show the range (upper figure) and standard deviation (lower figure) for the three seasonal periods for all data in each of the respective simulations. The values indicate the expected smoothing effect from pixel aggregation, though the variability is never totally removed. The Spring and Fall values demonstrate the highest variability. This variability is likely a function of the spatial heterogeneity seen in the irregular Landsat albedo frequency seen in Figures 2(b) and 2(d), and is expected since the Spring and Fall periods would present the most heterogeneous canopy as crops are planted, then either harvested or senesced. The Summer period has a smaller variability from a more fully developed canopy providing a more uniform signal. The uncertainty for all seasons is sufficient to indicate potential model influence at fine scales, up to at least 1000 m.

Figure 11. Albedo ranges (upper figure) and standard deviation (lower figure) for the Bondville site. The results for Bondville indicate a resolution of at least 1000 m to stabilize the effect from geometric distortion, though it is not fully removed even at 2000 m.

Several MODIS biophysical products are generated as discrete temporal composites, using only the single best quality per-pixel observation (e.g. cloud-free and/or minimum VZA) for a given composite period – usually 8–16 days. The series of plots in Figure 12 group the simulated data into discrete 8-day groups and calculate a range for each period. The 250 and 500 m groups demonstrate persistent variability greater than 0.01 albedo for all seasons, which is sufficient to influence an expected result within any given 8-day group, and suggest that care is needed when using observations from different geometries during a given composite period.

Figure 12. Discrete eight-day groups of simulated albedo for the Bondville site.

4.2. Harvard Forest LTER, MA

The variability in Figure 13 appears similar to the Bondville results for the 250 m resolution. However, the reduction in variability is achieved more rapidly than for Bondville, the majority occurring at 500 m, and a relative minimum reduction is achieved at 1000 m, similar to that for the Bondville site. The level of spatial heterogeneity at the Harvard site appears constrained to finer resolutions, and becomes more homogeneous at coarser resolutions. The Landsat albedo frequencies for the site (Figure 3) support this, particularly the Summer results that demonstrated a skewed distribution that is seen here as the highest level of uncertainty. The standard deviation values are small, suggesting the data are stable when using a mean value. However, the uncertainty remains a potentially significant contribution to the overall accuracy of the signal at finer resolutions.

Figure 13. Range (upper plot) and standard deviation (lower plot) for the Harvard Forest site. The albedo data for Harvard Forest demonstrate a continuing reduction of error as pixel size increases.

Figure 14 groups the simulated data into similar discrete 8-day groups as Figure 12. The Spring and Fall scenes demonstrate low variability from period to period and remain below for all resolutions, contrary to the global variability in Figure 13, which is cumulative over the range of all values in the simulation. The discrete results suggest the site is far more homogeneous and stable through geometry, though some care should be taken with fine-resolution data that present the highest per-group uncertainty. The results for Summer, however, demonstrate high levels of variation in the 250 m data. The effect is mitigated through aggregation, becoming more stable as the other seasons. The results for Summer suggest that 250 m data are sensitive to the GSD-to-grid assignment in developed tree canopies.

Figure 14. Discrete 8-day groups of simulated albedo for the Harvard site.

4.3. Konza Prairie LTER, KS

The results for Konza Prairie (Figure 15) demonstrate a low, and stable overall uncertainty for all resolutions. The standard deviation is insignificant, but the range data, while small, remain a constant source of significant uncertainty despite aggregation. The land cover seen in Figure 4(a) appears relatively uniform compared to Bondville and Harvard Forest, indicating a high level of spatial homogeneity for all seasons, and is supported in the Landsat albedo frequencies in Figure 4.

Figure 15. Albedo ranges and standard deviation for the Konza Prairie site. The relatively homogeneous land cover at the Konza Prairie site demonstrates an overall small and consistent uncertainty for all resolutions.

The discrete analysis (Figure 16) shows the same basic pattern of low uncertainty between groups similar to the Harvard site. The Fall data show the largest variation, though this is small. Overall, the results for Konza Prairie support MODIS data present a consistent uncertainty for homogeneous landscapes.

Figure 16. Discrete 8-day groups of simulated albedo for the Konza site.

4.4. Lake Tahoe, NV, and Jornada, NM

The final two sites are presented together since only one season was analyzed for each. The Lake Tahoe site was chosen for the interesting feature of having the exact site location defined in the lake, but the aggregation of pixels ultimately incorporating land features at coarser resolutions. Figure 5(a) illustrates the composition of pixels at increasing resolutions and the level of land features that are added to the signal. The level of variability is inverse to the expected smoothing for resolutions greater than 500 m seen in Figure 17, as the spatial variability of land features is introduced. Ultimately, the variability is the highest level for all the measured sites, but at the coarsest resolution. The integration of potentially disparate land covers into any given GSD is a consequence of the MODIS instrument's wide field of view. The variable land cover integration is carried forward into the gridded product, so care should be taken when making assumptions about per-pixel stability.

Figure 17. Albedo ranges and standard deviation for the Lake Tahoe and Jornada sites. The actual site center for Lake Tahoe is located offshore, but inclusion of signal from land features is steadily introduced at increasing spatial resolutions increasing the variability in the signal. The Jornada data present a stable level of uncertainty for all resolutions.

The Lake Tahoe site discrete results (Figure 18(a)) are similar to the global, in that nominal 250 m or aggregating data to 500 m is the most representative for this type of extreme heterogeneous land cover. Aggregation to coarser resolutions introduces a level of uncertainty to each period, unlike the other land covers represented in this research.

Figure 18. Discrete 8-day groups of simulated albedo for the Lake Tahoe and Jornada sites.

The Jornada site is relatively homogeneous, which can be seen in Figure 5 and is confirmed by the Landsat albedo frequencies in the figure. The homogeneity of the site is further confirmed with a low level of uncertainty for all resolutions in Figure 17.

The discrete analysis of the Jornada data (Figure 18(b)) shows that the 250 m data return the lowest uncertainty of all the resolutions. This result, and a similar finding for the Konza Prairie site, suggests that aggregation introduces heterogeneity from adjacent areas not otherwise incorporated in the nominal pixel. However, aggregation to resolutions needed to stabilize variability in other sites does not necessarily introduce a significantly different uncertainty.

4.5. Average resampling

Current MODIS land data products use the highest proportion of a GSD intersecting a given grid space, obscov, to assign a grid value (Wolfe, Roy, and Vermote 1998; Yang and Wolfe 2001). However, the resampling of a widely variable GSD to a fixed 250 m grid cell does not necessarily represent the sample space. Specifically, the misalignment of observations occurs when two observations partially intersect a grid cell and only a single value, with the largest coverage, is assigned; the feature of interest may be misassigned if the observation coverage is insufficient.

A second method of data processing was attempted as part of this study, to mitigate the inconsistent assignment resulting from the nearest neighbor (obscov in production) resampling. This method uses the average of all observations intersecting the grid cell to assign the grid a value. The average resampling technique is more representative of the changing GSD through the inclusion of all data informing a given grid space. Average resampling is processed from swath at all resolutions rather than aggregating the results of resampling a single resolution to progressively coarser resolutions as was done above. In a production environment, the simple average used here should be replaced with a weighted average using obscov, the same method used for aggregation in production. The results are seen in Figure 19(a).

Figure 19. Albedo profiles for all sites using an average resampling method for all resolutions (a). All sites except for Lake Tahoe demonstrate a reduction in range and standard deviations at fine resolutions when compared to the nearest neighbor resampling (b). The average resampling method presents a more consistent and smaller error at a finer resolution.

The immediate effect of the average resampling method is a reduction of range by at least half for the worst sites (from to ) and a reduction for other sites to generally albedo across all resolutions. The Lake Tahoe site remains an exception, demonstrating an increase in range beginning at a finer resolution. The standard deviation is similarly reduced for finer resolutions. However, some level of geometric error remains. The Bondville site for the Fall period demonstrates a clear increase at coarser resolutions, which is likely reflective of the introduction of spatially heterogeneous land cover similar to the Lake Tahoe site and the inherent sensitivity of averages to outlier values. The Lake Tahoe site does not benefit, and in fact suffers from the average resampling method.

These results indicate a more uniform choice of aggregating data to 1000 m to mitigate geometric effects and produce a more consistent error across land covers. The average resampling swath method also presents a better tradeoff between resolution and accuracy.

5. Discussion and conclusion

The objective of this research was to investigate the effect of geometric distortions inherent in the scanning MODIS system using broadband albedo as an example to quantify uncertainty. It should be noted that the purpose of this paper is not to assess the accuracy or model methods used for the MODIS standard albedo product (MCD43) or any other MODIS standard product. Rather, the findings of this work are meant to inform users of MODIS land products of the level of uncertainty added by the distortion inherent in MODIS data due to off-nadir viewing, using albedo as a meaningful physical variable. A quantified error budget associated with the spatial representativeness of wide field-of-view data such as MODIS will improve the accuracy and inform end use of current and future product generations through a better understanding of specific sources of error. The work has led to several valuable conclusions.

First, the assignment of a variable daily geometric footprint to a consistent geolocation grid has a significant influence on daily observations of synthetic albedo. Variability in albedo attributable to geometric effects alone ranged as high as –. Even applying a more robust assumption of maximum error that is half the range of the results still presents uncertainty in albedo for most sites at nominal resolution. These ranges from geometric and geolocation effects alone would be sufficient to influence climate models utilizing MODIS data. The standard deviation is much smaller for the same data and range from –, suggesting algorithms integrating multiple GSDs into a single, adjusted value suffer less from the inconsistent geometry.

Second, the research illustrates that error from geometric distortions can be mitigated through the smoothing of data by aggregation of pixels to coarser resolutions. Pixel resolutions in excess of 1000 m are required to mitigate the geometric influence from some land cover types, though 1000 m provides a relatively consistent result for all sites. The level of aggregation is consistent with previous findings using gridded data (Tan et al. 2006), as well as the known 4.8 times GSD growth along-scan.

The fact that complete removal of noise is not fully achieved even at coarse resolutions up to 8 times the nominal 250 m resolution when the known GSD of the instrument grows up to 4.8 times nominal resolution indicates that models using these data must incorporate an uncertainty attributable to geometric distortions before model error can be assessed. It should be noted that the current MCD43 BRDF/Albedo is generated at 500 m grid resolution rather than the 250 m geometries used for the above simulation (Schaaf et al. 2014).

The applicability of the results presented here is significant despite this difference in resolution for several reasons. The Red and NIR bands in the MCD43 product are nominal 250 m MODIS bands that are aggregated to 500 m using a more robust method than that presented here. The results for 250 m bands presented here clearly demonstrate that the smoothing from just a simple aggregation method to 500 m helps stabilize the geometric uncertainty, and the 500 m nominal bands could likewise benefit from aggregation to 1000 m. Additionally, the geometries for all MODIS detector resolutions function the same, but at different scales. Therefore, even the simple assumption of doubling of the minimum effective GSD for 250 m grid data (Campagnolo and Montaño 2014) at nadir (e.g. ) suggests that a 500 m grid product is not representative of the space being observed through varying geometry. While MODIS products are designed to accurately present the portion of the space observed in a given grid cell, there is no specific benefit from using a smaller grid cell size than a cell size more representative of the maximum variable observation footprint, for example; a resolution increasing the average obscov to greater than 30% (Tan et al. 2006) would be beneficial to stabilize the spatial signal.

Finally, the assessment of an alternate resampling method demonstrates averaging all swath pixel that fall into a grid space rather than nearest neighbor (as a proxy for obscov) is a more consistent spatial representation across the range of GSD. Future research of alternate data-processing methods should also investigate the efficacy of other equal area projections that can improve the swath-to-grid alignment over the sinusoidal projection. Campagnolo and Montaño (2014) demonstrate that the sinusoidal grid causes a degradation of the resolution for locations away from the Equator and the central meridian of the projection from inconsistent swath and grid orientation, which might be improved with a more localized projection.

The results demonstrate that the gridded data from different GSDs are not inherently interchangeable, and selection of data for any given period of geometry introduces a level of uncertainty that should be factored into product accuracy assessment for these data. This is particularly true not only of nominal 250 m data, but also for other resolutions. Spatial homogeneity for the site is an important factor when considering observation stability across the range of GSDs. Relatively homogeneous sites such as Konza Prairie, Jornada, and Harvard Forest provide a stable, though potentially significant, level of uncertainty for all resolutions.

Spatially heterogeneous sites such as Bondville can be smoothed through aggregation to similar levels of uncertainty to accommodate the spatial variability. Field size is also a factor when considering the most appropriate scale for MODIS data. Sites such as Bondville, with agricultural fields smaller than a given GSD, will present different portions of those fields with every orbit at finer resolutions. Thus, these land cover types benefit from a coarser scale to stabilize the signal. Conversely, some care should be taken not to aggregate the data to a level that no longer represents the land cover of interest, as the Lake Tahoe results demonstrated with two disparate adjacent land cover types.

The MODIS instruments have proven to be valuable contributions to the ongoing monitoring of Earth's surface (Gitelson, Kaufman, and Merzlyak 1996; Huete et al. 1999; Ji et al. 2008). There remains an important tradeoff between achieving high-temporal frequency observations and spatial accuracy; the results shown in this work clearly demonstrate the importance of minimizing measurement error arising from known instrument artifacts. However, the presence of known artifacts does not negate the use of existing nominal-scale data. Rather, the final use informs the level of utility.

Nominal 250 m gridded data are adequate when the final accuracy needed is relatively low as might be the case for a simple forest/non-forest classification or visual analysis. The results using standard deviation also indicate that the uncertainty of measurements is stable when considering products that integrate multiple GSDs into a single measurement. Physical variables or analysis assuming a stable spatial signal through time, per-pixel comparison from two different 8-day composites, for example, and requiring a high level of accuracy or measurement sensitivity can benefit from adding a geometric uncertainty factor based on the known variability of the GSD to increase the precision.

The use of whisk-broom sensors will continue and there is considerable interest in developing long-term consistent climate data records from satellites (Justice et al. 2011). The Visible Infrared Imaging Radiometer Suite (VIIRS) instrument, part of the orbiting Suomi NPP mission and to be flown on the upcoming Joint Polar Satellite System mission, will provide data that are aggregated along-scan to help control for the known effects of GSD growth. The intention of the aggregation is to provide a more consistent observation footprint than MODIS does currently, and early spatial performance (Wolfe et al. 2012) indicates the method is effective. Standard products from VIIRS, similar to the current MODIS suite, are available in beta form now. Future research will investigate efficacy of the cross-scan pixel aggregation using the methods presented here.

Acknowledgements

The author would like to thank Chris Justice and Eric Vermote, from the Department of Geographical Sciences at the University of Maryland, and Manuel Campagnolo from the Lisbon University, Institute of Agronomy, for their valuable comments and suggestions. Contributions to this work also came from many discussions with other faculty at UMD and NASA, in particular Mashihiro Nishihama at the Goddard Space Flight Center, whose expertise in simulation of MODIS data made this work possible. Finally, the author would like to thank the anonymous reviewers for their insight and comments.

Disclosure statement

No potential conflict of interest was reported by the authors.