Estimation of incident photosynthetically active radiation from
Moderate Resolution Imaging Spectrometer data
Shunlin Liang,
1
Tao Zheng,
1
Ronggao Liu,
2
Hongliang Fang,
1
Si-Chee Tsay,
3
and Steven Running
4
Received 3 October 2005; revised 14 March 2006; accepted 28 April 2006; published 8 August 2006.
[1] Incident photosynthetically active radiation (PAR) is a key variable needed by almost
all terrestrial ecosystem models. Unfortunately, the current incident PAR products
estimated from remotely sensed data at spatial and temporal resolutions are not sufficient
for carbon cycle modeling and various applications. In this study, the authors develop a
new method based on the look-up table approach for estimating instantaneous
incident PAR from the polar-orbiting Moderate Resolution Imaging Spectrometer
(MODIS) data. Since the top-of-atmosphere (TOA) radiance depends on both surface
reflectance and atmospheric properties that largely determine the incident PAR, our first
step is to estimate surface reflectance. The approach assumes known aerosol properties
for the observations with minimum blue reflectance from a temporal window of each
pixel. Their inverted surface reflectance is then interpolated to determine the surface
reflectance of other observations. The second step is to calculate PAR by matching the
computed TOA reflectance from the look-up table with the TOA values of the satellite
observations. Both the direct and diffuse PAR components, as well as the total shortwave
radiation, are determined in exactly the same fashion. The calculation of a daily average
PAR value from one or two instantaneous PAR values is also explored. Ground
measurements from seven FLUXNET sites are used for validating the algorithm. The
results indicate that this approach can produce reasonable PAR product at 1 km resolution
and is suitable for global applications, although more quantitative validation activities are
still needed.
Citation: Liang, S., T. Zheng, R. Liu, H. Fang, S.-C. Tsay, and S. Running (2006), Estimation of incident photosynthetically active
radiation from Moderate Resolution Imaging Spectrometer data, J. Geophys. Res., 111, D15208, doi:10.1029/2005JD006730.
1. Introduction
[2] Terrestrial ecosystems are thought to be a major sink
for carbon at the present time [Houghton et al., 1998].
Different terrestrial ecosystem models, including models of
terrestrial biogeochemistry [Running and Gower, 1991],
global vegetation biogeography [Prentice et al., 1992;
Prince, 1991; Prince and Goward, 1995], dynamic vegeta-
tion [Foley et al., 1996] and land-atmosphere exchange
processes [Bonan et al., 2002; Dai et al., 2003; Dickinson,
1995; Sellers et al., 1996], have been developed for the
function and dynamic nature of ecosystems, along with their
roles in the global carbon, nutrition, and water cycles.
Almost all these models contain the physiological processes
involved in photosynthesis and stomatal regulation that
control the exchange of water vapor and carbon dioxide
(CO
2
) between vegetation canopies and the atmosphere. To
initiate, calibrate and validate these models, various input
data sets from various sources, particularly from remote
sensing, are urgently required.
[3] Incident photosynthetically active radiation (PAR) is
one such key input and is needed to model photosynthesis
from single plant leaves to complex plant communities. The
net primary production or the rate of carbon accumulation
by terrestrial plant communities per unit area (C) is the
difference of gross photosynthesis by the canopy (A
g
) and
autotrophic respiration of the stand (R):
C ? A
g
C0 R ?1?
To use remote sensing products as key inputs, the
production efficiency models determine A
g
at varied spatial
scales [Field et al., 1995; Potter et al., 1993; Prince and
Goward, 1995; Running et al., 1999, 2000]:
A
g
/
Z
et??C1fPAR t??C1PAR t??dt ?2?
where e(t) is the radiation efficiency factor, fPAR represents
the fraction of the PAR absorbed by green vegetation. A
JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 111, D15208, doi:10.1029/2005JD006730, 2006
Click
Here
for
Full
Article
1
Department of Geography, University of Maryland, College Park,
Maryland, USA.
2
Institute for Geographical Science and Natural Resources Research,
Chinese Academy of Sciences, Beijing, China.
3
Climate and Radiation Branch, NASA/Goddard Space Flight Center,
Greenbelt, Maryland, USA.
4
School of Forestry, University of Montana, Missoula, Montana, USA.
Copyright 2006 by the American Geophysical Union.
0148-0227/06/2005JD006730$09.00
D15208 1of13
number of process-based global NPP models do not
explicitly calculate fPAR, but compute a leaf area index
(LAI) [Ruimy et al., 1999]. LAI and fPAR have a simple
functional relationship. Monteith [Monteith, 1972, 1977]
considered R to be a constant fraction of A
g
, and, in such a
case, spatial and temporal variations of C will be strongly
determined by spatial and temporal variations of the
absorbed PAR (APAR)(APAR = fPAR C1 PAR). fPAR is being
produced by the MODIS (Moderate-Resolution Imaging
Spectroradiometer) team as an Earth Observing System
(EOS) standard product [Myneni et al., 2002], but high-
resolution incident PAR product over land simply does not
exist although some efforts are being made [Van Laake and
Sanchez-Azofeifa, 2004, 2005].
[4] A worldwide observation network for measuring
incident PAR has not yet been established. The only
practical means of obtaining incident PAR at spatial and
temporal resolutions appropriate for most modeling appli-
cations is through remote sensing. Frouin and Pinker
[Frouin and Pinker, 1995] reviewed the methods for esti-
mating incident PAR from International Satellite Cloud
Climatology Project (ISCCP) and Total Ozone Mapping
Spectrometer (TOMS) observations. There have been two
monthly ISCCP PAR products at 2.5C176 resolution. The
ISCCP-BR PAR product was calculated from incident total
shortwave flux divided by 2 [Potter et al., 1993] using the
ISCCP C1 data that are composed of both geostationary and
polar-orbiter observations. The ISCCP-PL PAR product was
derived from the PAR algorithm of Pinker and Laszlo
[Pinker and Laszlo, 1992] also using the ISCCP C1 data.
The TOMS PAR product was based on the method of Eck
and Dye [Eck and Dye, 1991] using ultraviolet reflectivity.
Dye and Shibasaki [1995] compared these three PAR
products and found that ISCCP-BR and ISCCP-PL are
typically 12?16% and 8?12% greater than TOMS PAR,
respectively. Comparison of these three PAR products with
ground data at a midlatitude site for nonwinter months by
the same authors yielded the root mean square (RMS)
differences of 28.1%, 13.7%, and 7.2%, and biases of
+25.9%, +12.0%, and +2.8%, respectively.
[5] Because the high-resolution incident PAR over land
surface is not a standard EOS product, the MODIS science
team has to disaggregate the NASA data assimilation PAR
product of 3-hourly 1C176 by 1.5C176 spatial resolution to produce
1 km NPP and net photosynthesis (PSN) products using the
BIOME-BGC model [Running et al., 1999, 2004]. The
GLObal Production Efficiency Model (GLO-PEM) has to
use 1C176 by 1C176 monthly TOMS PAR product for calculating
global 8 km NPP product [Prince and Goward, 1995]. The
TOMS PAR product is also limited to 66C176Nto66C176S latitude
and is not available after 1989. The simplified PAR flux
models are being used to produce the incident PAR flux
over ocean from MODIS data [Carder et al., 2003] and Sea-
viewing Wide Field-of-View Sensor (SeaWiFS) data
[Frouin et al., 2000]. No PAR product has been routinely
generated from MODIS or SeaWiFS over land. The incident
PAR product constitutes a major uncertainty in carbon cycle
modeling, and further developments are urgently needed.
[6] Satellite observations acquired at the top of the
atmosphere (TOA) contain information of both atmosphere
and surface. The current methods [e.g., Carder et al., 2003;
Frouin et al., 2000; Gu et al., 2004; Pinker and Laszlo,
1992; Pinker et al., 2003] assume that either the atmo-
spheric information is available from other sources (e.g., the
ISCCP PAR product with atmospheric climatology data as
input) or that water surface reflectances are known (e.g.,
SeaWiFS and MODIS PAR products over ocean).
[7] If the atmospheric parameters and surface spectral
albedos are known, the existing algorithms [e.g., Van Laake
and Sanchez-Azofeifa, 2004, 2005] may be directly applied.
However, these models are based on the simple two-stream
approximations for multiple scattering or even simpler
schemes. Computation is speeded up significantly with the
simplified models, but they usually are not accurate in
calculating multiple scattering that dominates when the
atmospheric optical depth is large (hazy or cloudy) and/or
the surface is very bright (snow/ice). Simple models are
easy to implement, and suitable for cases where no accurate
input parameters are available (e.g., processing GOES data
alone). Increasing the model complexity will increase dif-
ficulties of implementation for the regional and global
applications. Instead, the look-up table approach is used
in this study. It is specifically designed to achieve the high
accuracy and keep a simple form of the model. The look-up
table approach has been widely used in various scientific
investigations.
[8] Another issue is the separation of direct and diffuse
PAR radiation. The volume of shade within vegetation
canopies is reduced by more than an order of magnitude
on cloudy and/or very hazy days compared to clear, sunny
days because of an increase in the diffuse fraction of the
solar radiance. In a recent study, Gu et al. [2002] found that
(1) diffuse radiation results in higher light use efficiencies
by plant canopies; (2) diffuse radiation has much less
tendency to cause canopy photosynthetic saturation; (3) the
advantages of diffuse radiation over direct radiation increase
with radiation level; and (4) temperature as well as vapor
pressure deficit can cause different responses in diffuse and
direct canopy photosynthesis, indicating that their impacts
on terrestrial ecosystem carbon assimilation may depend on
radiation regimes and thus sky conditions. These findings
call for different treatments of diffuse and direct radiation in
models of global primary production, and studies of the
roles of clouds and aerosols in global carbon cycle [Gu et
al., 2003]. In fact, many land surface process models, such
as SiB2 [Sellers et al., 1996] and CLM [Dai et al., 2003],
have separated the direct and diffuse solar radiation. How-
ever, none of the existing PAR products separate these two
components.
[9] In this study, the authors develop a new method for
estimating incident PAR from MODIS data by estimating
both surface and atmospheric information simultaneously.
The key concept is the use of multitemporal signatures of
MODIS data. Since the basic requirement of this method is
the high temporal resolution, it potentially can be used for
many different polar-orbiting sensors, such as Medium-
Resolution Imaging Spectrometer (MERIS), SeaWIFS,
Advanced Very High Resolution Radiometer (AVHRR),
Global Land Imager (GLI) and Visible/Infrared Imager/
Radiometer Suite (VIIRS), and geostationary sensors, such
as geostationary operational environmental satellites
(GOES).
[10] MODIS is one of the sensors in NASA?s EOS Terra
platform launched on 18 December 1999 for imaging in the
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
2of13
D15208
morning, and the Aqua platform launched on 4 May 2002
for imaging on the afternoon. MODIS has 36 spectral bands
spanning from the visible to thermal-Infrared (IR) spectrum.
The spatial resolution varies from 250 m (bands 1 and 2),
500 m (bands 3?7), to 1000 m (bands 8?36). The swath
width of MODIS is about 2300 km, with the across-track
field-of-view angle of 110C176. Each MODIS sensor images the
earth on a 2-day repeat cycle, with a 1-day or more frequent
repeat at higher latitudes greater than 30C176 because of orbital
convergence. The details of the sensor characterization are
available elsewhere [Salomonson et al., 1989].
[11] The details of the new method are presented in
section 2. The validation results are shown in section 3. A
series of related issues are discussed in section 4. A brief
summary is given at the end of this paper.
2. New Method Description
[12] Incident PAR depends mainly on the atmospheric
properties, but also to a lesser extent on surface reflectance.
It can be demonstrated by the following formula for
calculating downward spectral flux F(m
0
) at the solar
zenith angle q
0
(m
0
= cos(q
0
)) over a Lambertian surface
[Chandrasekhar, 1960; Liang, 2004]:
F m
0
???F
0
m
0
???
r
s
C22r
1 C0 r
s
C22r
m
0
E
0
g m
0
?? ?3?
where F
0
(q
0
) is the downward flux without any contribu-
tions from the surface, r
s
is surface reflectance, C22r is the
spherical albedo of the atmosphere, E
0
is the extraterrestrial
solar irradiance, g(m
0
) is the total transmittance (direct plus
diffuse) in the solar illumination direction. All radiometric
variables are the function of wavelength. PAR is the
integration of the downward spectral fluxes integrated from
400 nm to 700 nm:
PAR m
0
???
Z
700
400
F
l
m
0
??dl ?4?
[13] We use the energy unit Wm
C02
in our look-up tables.
Note that many ecosystem process models typically employ
PAR data in quantum units (photosynthetic photon flux
density, mmol m
C02
s
C02
) and the conversion between the
energy units to the quantum units was recently discussed
by Dye [2004]. It is clear that the surface contribution to the
incident PAR mainly depends on r
s
C22r. If the atmosphere is
very clear (i.e., C22r is small) or the surface reflectance r
s
is
low, the surface contribution to the incident PAR is rela-
tively small.
[14] Our new method first estimates surface reflectance
from multitemporal imagery and then PAR flux for each
imagery. The basic procedure is composed of two steps,
including (1) determination of the surface reflectance from
the ??clearest?? observations during a temporal window and
(2) calculation of incident PAR from the determined surface
reflectance and TOA radiance/reflectance using the table
look-up approach. Determining the length of the temporal
window needs to be done carefully. Obviously, it must be
short enough so the surface properties do not change
dramatically, but long enough to include adequate clear
observations. In all our case studies, it seems a period of one
to three months is a reasonable choice. The details in each
step are discussed below after introducing how to create
look-up tables.
[15] Note that this new method differs from other meth-
ods in that surface reflectance and atmospheric properties
are estimated simultaneously from imagery itself. For
MODIS imagery, both aerosol/cloud optical depth and
surface reflectance are the standard products from the
MODIS science team. We can use these products for
calculating PAR [e.g., Van Laake and Sanchez-Azofeifa,
2004, 2005]. However, some of these algorithms are not
well developed at this point. For example, the aerosol
optical depth products over land has many gaps (no retrieval)
and is only available when the surface is densely vegetated
[Kaufman et al., 1997a, 1997b]. Besides, the MODIS
atmospheric products have much coarser spatial resolution
(C2410 km). Therefore the method of Van Laake and
Sanchez-Azofeifa is not suitable for estimating incident
PAR from MODIS data at 1 km resolution.
2.1. Creating Look-Up Tables (LUT)
[16] ForaLambertiansurface,theTOAupwellingradiance
can be calculated from the classic formula [Chandrasekhar,
1960; Liang, 2004]:
I m
0
; m; f???I
0
m
0
; m; f???
r
s
1 C0 r
s
C22r
m
0
E
0
g C0m
0
??g m?? ?5?
where I(m
0
, m, f) is the observed TOA radiance at the
viewing zenith angle q, m = cos(q), and relative azimuth
angle f, I
0
(m
0
, m, f) is path radiance without the surface
contributions, and g(m) is the total transmittance from the
surface to the sensor.
[17] Instead of online calculation of these atmospheric
functions in equation (5) on the pixel basis that can be
computationally very expensive, they are tabled offline as a
function of aerosol or cloud optical depth at different
viewing geometries. MODTRAN [Berk et al., 1998] is used
in this study, since it has been explored extensively in
various simulations herein. In MODTRAN, multiple scat-
tering is accurately calculated by DISORT [Stamnes et al.,
1988]. For each aerosol model at each specific sun-viewing
geometry with a range of atmospheric visibility, three
surface reflectance values are specified (0.0, 0.5 and 0.8)
for the range of wavelengths at the interval of 15 wave
numbers. MODTRAN outputs at each wavelength interval
for these three surface reflectance are used to solve for
F
0
(m
0
), F
d
(m
0
)=m
0
E
0
g(m
0
), C22r, I
0
(m
0
, m, f), and g(m). These
Table 1. A Typical Look-Up Table for Each Aerosol Model at a
Specific Sun-Viewing Geometry, Where i = 1, 3, 4, v Denote
MODIS Three Bands: 1 (Red), 3 (Blue), 4 (Green), and the Whole
Visible Spectrum
a
Visibility
(Vi), km
I
0
(i)
i = 1, 3, 4
F
d
(i)
i = 1, 3, 4, v
g
(i)
(m)
i = 1, 3, 4
C22r
(i)
i = 1, 3, 4, v
C1C1C1 C1C1C1 C1C1C1 C1C1C1 C1C1C1
C1C1C1 C1C1C1 C1C1C1 C1C1C1 C1C1C1
a
Direct and diffuse PAR and total shortware radiations are three separate
columns in the look-up table.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
3of13
D15208
quantities are then integrated with the MODIS sensor
spectral response functions over wavelength to get the
average values for each wave band. The average values of
F
0
(m
0
), F
d
(m
0
) and C22r over the visible spectrum (PAR region)
also are calculated. Thus, for each aerosol model at a
specific sun-viewing geometry (m
0
, m, f), there is a table
that looks like Table 1. For clouds, the table looks the
similar. Instead of changing visibility, cloud extinction
values are changed.
[18] In this study, the following values have been used in
the MODTRAN simulations: aerosol types (rural, maritime,
urban and tropospheric), cloud types (altostratus, stratus,
stratus/stratocumulus, and nimbostratus), solar zenith angle
(0C176,20C176,40C176,50C176,60C176,65C176,70C176,75C176, and 80C176), viewing
zenith angle (0C176,15C176,30C176,45C176, and 65C176), relative azimuth
angle (0C176,30C176,60C176,90C176, 120C176, 150C176, and 180C176), and
visibility for characterizing aerosol loadings (5, 10, 20,
30, 50 and 100 km). The default values of cloud thickness
and base heights are used, but extinction coefficient varies
for five values specified at 550 nm for each cloud type (see
Table 2).
[19] The impacts of absorptive gases on incident PAR are
examined through MODTRAN simulations. Because of
their spatial and temporal variations, simulations are limited
only to ozone and water vapor. The integrated transmittance
of the PAR region at nadir and their concentrations are
shown in Tables 3 and 4. By ignoring their variations,
the uncertainty of the resulting PAR is clearly quite
small. Therefore they are not considered as variables in
the look-up tables.
[20] The look-up table approach is simple and fast, but
has its limitations. It works only when the output is fairly
linear with respect to changes in independent variables and
the interactive effects of different variables are weak. It also
limits the number of variables that can be considered in the
simulation. The intervals of each variable in the look-up
table need to be determined through a series of sensitivity
studies so that a linear interpolation in the look-up table
searching is valid.
2.2. Determination of Surface Reflectance
[21] The MODIS science team is routinely producing the
land surface reflectance product routinely. One of its major
inputs is the aerosol optical depth, which is estimated by
using the ??dark-object?? method. The ??dark-object?? method
has been widely used [Liang, 2004], but its limitation is that
it is suitable only for dense vegetation canopies. In this
study, we have developed a simpler method since incident
PAR is much less dependant on surface reflectance.
[22] The central idea of this method is to detect the
??clearest?? observations in the temporal dimension for each
pixel. Assuming that the aerosol optical depths for the
??clearest?? observation are known, its surface reflectance
can be determined by assuming the surface is a Lambertian
reflector. Surface reflectances of other ??hazy/cloudy??
observations are then interpolated from these ??clearest??
observations.
[23] The ??clearest?? observations are identified on the
basis of their minimal blue band values. The underlying
assumption is that when the atmosphere becomes more
turbid, the blue band TOA radiance will be larger. On the
other hand, surface reflectance at blue wavelengths is
usually the lowest in the visible region [Eck and Dye,
1991]. Thus the blue band is a good choice for us to
identify the ??clear?? atmospheric conditions. To avoid
picking up the shadowing observations, we convert the
TOA radiance of the blue band into the apparent surface
reflectance from the look-up table values with a known
aerosol visibility value for a very clear atmospheric condi-
tion (100 km is assumed in this study). If an observation
contains cloud shadows, the converted surface reflectance
value is usually negative and therefore flagged and excluded
in determining surface reflectance. Here a fixed percentage
of the total observations within the temporal window is used
to determine the number of clear observations (e.g., 10%) at
this stage. Ideally, it should be variable both spatially and
temporally. More experiments may be needed in the future.
[24] As soon as the surface reflectance of these ??clear??
observations is calculated, the surface reflectance of other
??hazy?? observations can be linearly interpolated from them.
To ensure that sufficient amount of ??clear?? observations are
included in the temporal window, the length of the temporal
window should be long enough. Thus surface reflectance
within the temporal window may change and an interpola-
tion process is necessary. In our experiments, it is found that
the length of the temporal window of 1?3 months is
reasonable.
Table 2. Cloud Extinction Coefficients (km
C01
) at 0.55 mm, Thickness and Base Heights
Altostratus
Cloud
Stratus
Cloud
Stratus/Stratocumulus
Cloud
Nimbostratus
Cloud
Extinction coefficients 1 1 1 1
Extinction coefficients 5 5 3 5
Extinction coefficients 20 15 10 10
Extinction coefficients 50 30 15 30
Extinction coefficients 128.1 56.9 38.7 45
Thickness, km 0.6 0.67 1.34 0.5
Base height, km 2.4 0.33 0.66 0.16
Table 3. Ozone Concentration and the Integrated Transmittance
Over the PAR Region
Concentration, gm/m
2
Transmittance
0.0 1.0
0.6393 0.9986
3.4804 0.9919
7.0319 0.9838
10.5833 0.9757
14.1347 0.9679
17.6862 0.9603
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
4of13
D15208
[25] Although PAR flux is less dependant on surface
reflectance, we need to consider an abrupt significant
change of surface reflectance. The extreme case is the
snowfall and snowmelt. The snow observations have to be
labeled and separated from nonsnow observations. They
must be processed separately but in the same fashion. The
MODIS team is routinely producing snow and ice cover
map. These maps can be used to distinguish snow and
cloud, which has been an issue in previous studies [e.g.,
Pinker et al., 2003]. What we do in this study is to group
snow and nonsnow observations using the MODIS snow
cover maps and apply the same procedures to these two
groups separately.
2.3. Searching LUT for Calculating Incident PAR
[26] Given the surface reflectance from section 2.2 and
TOA reflectance that is converted from digital numbers
using the calibration coefficients, we can search the look-up
tables and calculate incident PAR.
[27] For each pixel, if we know surface reflectance r
s
(i)
,
i = 1, 2, ... N, the TOA radiance
^
I
j
(i)
for these N bands
can be predicted for each raw j in Table 1 using equation (5).
The observed TOA radiance I
(i)
are then used to compare
with
^
I
j
(i)
for determining the visibility value using linear
interpolation since each j corresponds to one visibility
value. The visibility value is the intermediate quantity and
mainly used to interpolate other variables in the look-up
tables. For each band (i), we can determine a visibility
value. Because of the all possible uncertainties due to
calibration and various assumptions, the resulting visibil-
ity values from different visible bands are usually differ-
ent. Several schemes have been explored to calculate the
PAR value.
[28] In scheme 1, for each estimated visibility value from
each band, the downward flux of each band (F
i
) is calcu-
lated by equation (3) with the quantities interpolated from
Table 1 using the estimated visibility value. The PAR value
is the linear combination of the spectral downward fluxes:
PAR ? a
0
?
X
N
i?1
a
i
F
i
?6?
where a
i
are the coefficients determined from simulations.
[29] In scheme 2, there is an estimated visibility value for
each band, i, that is then used for interpolating the average
values of several quantities (F
0
(v)
, C22r
(v)
, F
d
(v)
) from Table 1 so
that PAR can be calculated in an aggregated fashion:
PAR
i??
m
0
???F
v??
0
m
0
???
r
v??
s
C22r
v??
1 C0 r
v??
s
C22r
v??
F
v??
d
m
0
?? ?7?
where the broadband visible reflectance r
s
(v)
is calculated
from the spectral reflectance on the basis of our previous
work [Liang, 2001, 2004]. For MODIS, the formula is
r
v??
s
? 0:331r
1??
s
? 0:424r
3??
s
? 0:246r
4
s
?8?
Thus the final PAR value is the average value over these
bands
C18
PAR =
1
N
P
N
i?1
PAR
(i)
C19
.
[30] Scheme 3 is similar to scheme 2 except that the
averaged visibility value is first calculated from all three
bands and then the PAR is calculated using equation (7).
[31] Note that the direct, diffuse and total PAR are three
columns in the look-up table, and each visibility value
corresponds to these three quantities. As long as the
visibility value is determined, they can be calculated in
the same fashion.
3. Data Analysis
[32] In the following examples, we will demonstrate how
the incident instantaneous PAR can be produced from
MODIS data using this look-up table approach. In the first
part, the validation results using ground measurements of
seven stations of FLUXNET will be presented. The MODIS
aerosol and cloud products are not used in these case
studies, since the estimation of incident PAR is not very
sensitive to the aerosol and cloud models. In the second
part, mapping incident PAR over the greater Washington, D.
C. region is presented to demonstrate that this method is
suitable for regional and global applications.
3.1. Validation
[33] The validation experiments were conducted
using ground measurements of incident PAR from seven
FLUXNETsites: Black Hill, Fort Peck, Lost Creek, Howard
forest (main tower), Oak Ridge, Sante Rem and Vaira
Ranch. For each site, a 3 km
*
3 km window (9 pixels)
of the MODIS TOA radiance (MOD02) and angular values
were extracted from the MODIS data sets ordered from the
EOS gateway. The ground measurements collected every
half-hour or 1 hour were compared with the retrieved
values. The measurement values closest to the time of
MODIS data acquisition were used to compare with the
value of the central pixel in the extracted 9-pixel window
without any interpolation. Because of the possible mismatch
in both space and time, it is difficult to characterize the
correlation between them using the traditional statistical
analysis. Another mismatch problem is the amount of
clouds and shadows existing in the different illumination
and viewing paths. This effect is illustrated in Figure 1,
showing the cloudy illumination path and clear viewing
path (Figure 1a), and the clear illumination path but cloudy
viewing path (Figure 1b). In the second case, the ground
measured PAR could be much larger than the retrieved PAR
from satellite observations even if the inversion is perfect.
The ideal solution is to employ three-dimensional atmo-
spheric radiative transfer modeling to account for all these
factors. However, incorporating such a sophisticated radia-
tive transfer calculation into the practical product generation
continues to be a challenging topic in remote sensing.
Table 4. Water Vapor Concentration and the Integrated Transmit-
tance Over the PAR Region
Concentration, gm/m
2
Transmittance
0.0 1.0
0.3703 0.9993
2.0162 0.9957
4.0735 0.9919
6.1309 0.9885
6.7258 0.9875
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
5of13
D15208
[34] Instead, the statistics are determined using robust
regression analysis technique, called the least trimmed
squares (LTS) regression [Rousseeuw, 1984].This regression
method minimizes the sum of the k smallest squared
residuals (n/2 < k < n) where n is the number of data
points. In this study, the largest squared residuals are
trimmed by 10% (i.e., k/n = 0.9). The fitted statistics of
all seven sites are summarized in Table 5a and shown in all
validation figures. For comparison, the statistics from the
ordinary regression analysis are summarized in Table 5b.
In Table 5b, we also provide the mean relative error,
defined as
1
N
P
N
i?1
1 C0
^y
i
y
i
C12
C12
C12
C12
C12
C12, where N is the total number of
measurements, y
i
and ^y
i
are the observed and estimated
PAR, respectively. The detailed analysis for each site is
given below.
Figure 1. Illustrations of (a) the cloudy illumination path
and clear viewing path and (b) the clear illumination path
but cloudy viewing path.
Table 5a. Summaries of Robust Regression Analysis at Seven
Validation Sites
Sites Intercept Slope
Scale of the
Residual
Robust
Multipled R
2
Black Hills C07.75 1.02 5.57 0.997
Fort Peck 3.57 1.00 36.33 0.996
Lost creek 2.39 1.00 19.15 0.999
Howard Forest 49.48 0.78 91.23 0.959
Oak Ridge C00.03 1.00 0.79 1.000
Santa Rem 10.05 0.99 49.79 0.976
Vaira Ranch C01.21 1.00 14.74 0.999
Table 5b. Summaries of Ordinary Regression Analysis at Seven
Validation Sites
Sites Intercept Slope
Residual
Standard
Error Multipled R
2
Average
Relative
Error, %
Black Hills 128.60 0.8645 116.6 0.9013 6.61
Fort Peck 193.13 0.8385 191.5 0.8804 15.1
Lost creek 66.801 0.9227 167.5 0.9247 10.2
Howard Forest 99.60 0.7394 161.8 0.8940 21.9
Oak Ridge C03.513 1.0202 15.22 0.9881 4.1
Santa Rem 182.59 0.8195 136.9 0.8661 7.6
Vaira Ranch 12.523 0.9773 85.65 0.9788 4.2
Figure 2. Comparisons of the retrieved PAR from three
visible bands: (a) green and blue and (b) red and blue.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
6of13
D15208
[35] The first site is Black Hill in North Dakota, which
has an evergreen needleleaf forest and temperate landscape.
The tower is located at 44C1769
0
28.8
00
N, 103C17639
0
W. The half-
hour measurement data in 2002 from days 182?304 were
used. Three schemes in searching the look-up tables are
mentioned in section 1.3. Different schemes produce slightly
different retrievals, but the differences are quite small. We
prefer to scheme 2. Further data analysis indicate that using
even only one band is sufficient to produce the PAR, instead
of the average estimates from all three bands, since the
retrievals from all three bands are very similar (see Figure 2).
These assessments can reduce the computation time and
space requirement significantly. Finally, the blue band
(band 3) is selected in this study.
[36] It was also found that selecting different cloud types
and aerosol types does not significantly affect the PAR
retrieval. Figure 3 shows the results of three different
aerosol types and two cloud types. Determining the aerosol
and cloud types from MODIS data is still very challenging.
Given the relative small differences, rural aerosol and
stratus cloud are used in the rest of the paper.
[37] The second site is the Fort Peck (48C17618.473
0
N and
105C1766.032
0
W), Montana station, located on the Fort Peck
Tribes Reservation, approximately fifteen miles north of
Figure 3. Comparisons of the retrieved incident PAR with three different aerosol types and two cloud
types.
Figure 4. Validation results at FLUXNET site Fort Peck.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
7of13
D15208
Poplar, Montana. A tower of 3 m high was installed in
November of 1999 over grassland with a surface elevation
of 634 m. The comparison of the retrieved and measured
half-hour incident PAR in 2002 (the whole year) is shown in
Figure 4. The overall fitting is very good, although there are
some pixels of overestimation with high PAR values and
underestimation with low PAR values.
[38] Note that all measurements are in the unit of
mmol/m
2
/s, but our look-up table PAR outputs are in Wm
-2
.
Forquickcomparisons,wesimplyusedtheconversionfactor
4.6.However,theconversionfactorissupposedtovaryunder
different conditions (see Dye [2004] for details).
[39] The third site is Lost Creek, (46C1764.9
0
N, 89C17658.7
0
W),
located in northern Wisconsin. The elevation is about 480 m
above sea level. The biome is mixed forest with canopy
height of 1?2 m. This site was established on 4 September
2000. The measurement data from both 2000 and 2001 were
used here and the validation results are shown in Figure 5.
The fitting is better than the two previous sites.
[40] The fourth site is Howland forest (Main Tower),
located at Howland, Maine (45.204C176N, 68.74C176W). Topo-
graphically, the region varies from flat to gently rolling with
a maximum elevation change of less than 68 m, The
landscape is composed of a deciduous evergreen needle
forest, boreal/northern hardwood, old coniferous, hemlock,
douglas fir, and evergeen coniferous. It is chiefly cold,
temperate and humid. The region is covered by the snow
of up to 2 m from December through March. Measurement
data in 2002 from dates 74?365 were downloaded. The
validation results are shown in Figure 6, which are the worst
in all cases since the retrievals are biased. We had difficulty
in figuring out what causes this bias, which is not obvious in
other sites. The instrument calibration could be a factor, but
we are making efforts to intercompare the PAR products
from other sensors.
[41] The fifth site is the Tapajos National Forest (Santa
Rem ? Km83) over the logged evergreen tropical broadleaf
forest at Brazil. The geographic coordinates of the tower is
3C1761
0
4.9
00
S, 54C17658
0
17.166
00
W. The measured PAR data started
from 1999, but only data in 2001 were used in this study.
The comparison is shown in Figure 7. The overall validation
results are good, although there exist some pixels of
Figure 5. Validation results at FLUXNET site Lost Creek.
Figure 6. Validation results at FLUXNET site Main
Tower.
Figure 7. Validation results at FLUXNET site Santa Rem.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
8of13
D15208
overestimation when the PAR values are large. Since this
site is in the tropic region, broken clouds and rapid move-
ments of cloud could be a contributing factor.
[42] The sixth site is the Walker Branch Watershed
(35C17657
0
31.56
00
N, 84C17617
0
14.76
00
W ), located at Oak Ridge,
Tennessee. It has a forest stand over 50 years old, having
regenerated from agricultural land. It has a temperate
climate, with a mixed species, broad-leaved forest, decidu-
ous forest, oak and hickory. This site is the only one where
both direct and diffuse radiation values are measured. The
validation results are shown in Figure 8. The direct com-
ponent has a much larger scatter, including a group of
underestimation pixels, but the agreements for the diffuse
and the total PAR appear much better. Because of the spatial
Figure 8. Validation results at FLUXNET site Oak Ridge:
(a) diffuse, (b) direct, and (c) total PAR.
Figure 8. (continued)
Figure 8. (continued)
Figure 9. Validation results at FLUXNET site Vaira
Ranch.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
9of13
D15208
and temporal mismatch between the ground measurements
and MODIS data, the large difference of the direct compo-
nent is not surprising. In the future validation, we ought to
investigate whether the inversion method also contributes to
the differences. For example, the direct PAR may be very
sensitive to the existence of partial clouds, while our
method assumes plane-parallel atmospheric properties. An-
other factor could be the conversion constant that coverts
the energy unit in our code to the quantum unit in the
measurement. We used a constant, but it could be variable
[Dye, 2004].
[43] The last site is the Vaira Ranch (38C17624.4
0
Nand
120C17657.044
0
W) at Ione, California, situated in a topograph-
ically undulating oak and grass savanna biome in eastern
California at the foot of the Sierra Nevada Mountains. The
surrounding landscape is grazed grass rangeland. The ele-
vation at the tower location is 129 m. A year of half-hour
measurements in 2002 was used in this study. The fitted
results are shown in Figure 9, and the agreement is
excellent.
[44] It is clear that this method can produce reasonably
accurate incident PAR products from MODIS observations.
However, there are large scatters in some cases. The spatial
and temporal mismatches between the ground ??point??
measurements and the MODIS pixel values are probably
the major source of these differences. Other factors
Figure 10. True color composite MODIS imagery of greater Washington, D. C., area in 8 days (22 May,
25 May, 29 May, 31 May, 5 June, 7 June, 8 June, and 10 June 2003).
Figure 11. Retrieved incident PAR maps from the corresponding MODIS images in Figure 10.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
10 of 13
D15208
may include the plane-parallel atmospheric radiative trans-
fer model that is used for creating the look-up tables
but fails to account for broken clouds in the atmosphere,
surface Lambertian assumption, linear interpolation in table
searching, uncertainty in determining surface reflectance,
and so on.
3.2. Mapping INCIDENT PAR Over the Greater
Washington, D. C. Area
[45] To test the method for the regional mapping, the
incident PAR over the greater Washington, DC area was
mapped in 2003. It is impractical to show all these images
here and the estimated PAR maps; instead, the images of the
following 8 days (142, 145, 149, 151, 156, 158, 159, and
161) are presented in this paper with variable atmospheric
conditions and visual spatial patterns. The true color com-
posite images (using bands 1,3 and 4) are shown in
Figure 10, and the estimated corresponding incident PAR
maps are shown in Figure 11. The image dimension is 800
by 800 pixels. It is clear there are very good correlations
between the original image, and the mapped PAR and the
patterns match very well.
4. Calculating Daily PAR
[46] Instantaneous PAR can be very useful to NPP and
other ecosystem models. However, many models have the
daily time step and thus the daily PAR product is more
desirable. If the high-frequency data from a geostationary
satellite (e.g., GOES) are available, calculating daily values
would be easier. An alternative method is explored in this
study.
[47] The daily PAR is calculated from two instantaneous
PAR products (MODIS/Terra (AM) and MODIS/Aqua
(PM). Experiments show that the diurnal cycle can be
approximated by cosine function of the local solar zenith
angle [e.g., Tarpley, 1979]. Following the same logic, the
morning and afternoon diurnal cycles are expressed by two
different functions. Specifically, MODIS/Terra instanta-
neous PAR values are used for determining the coefficient
of the morning function, and MODIS/Aqua instantaneous
PAR values are used for determining the afternoon function.
The daily PAR is simply the integration of these two
functions. The simpler version of this method is to conduct
the regression between MODIS-estimated instantaneous
values and the daily average.
[48] After analyzing these ground measurement data
(section 3.1), it was found the predicted daily mean PAR
flux, using one or two instantaneous MODIS PAR fluxes,
can be reasonably accurate. Figure 12 shows the fitting
results using in situ data of these FLUXNET sites. The
corresponding equations for calculating the daily 24-hour
average PAR from one morning PAR at 1030 local time
(LT) (PAR
am
) or one afternoon PAR (PAR
pm
) at 1330 LT or
both:
PAR daily ? 39:86 ? 0:0784*PAR
am
? 1:5277*10
C04
*PAR
2
am
PAR daily ?C06:637 ? 0:3571*PAR
pm
PAR daily ? 29:7 ? 0:0129*PAR
am
? 1:1715*10
C04
*PAR
2
am
? 0:128PAR
pm
8
>
>
<
>
>
:
There is a high degree of variation. Note that in the high
latitude regions there are multiple MODIS observations
each day from which the daily average values will be
predicted more accurately.
5. Incident PAR and Shortwave Insolation
[49] Since direct measurement of incident PAR is just
beginning and is still quite limited in space, it has been
assumed in various applications that incident PAR is half of
the incident shortwave radiation (insolation) that has rou-
tinely been measured by weather stations. Several studies
have indicated that this ratio is not a constant. Jacovides et
al. [2003] found that this ratio varies from 0.460 to 0.501
from the hourly measured values.
[50] Alados et al. [1996] found that this ratio also has
seasonal and daily variations. The ratio is smaller and more
variable in winter than in summer, and is usually low at
noon.
[51] On the other hand, the total shortwave radiation is
also required by any land surface process models and many
other applications. Using our new method, the total short-
wave downward flux can be organized in the same look-up
tables as incident PAR. Thus they can be determined
simultaneously. The surface reflectance of broadband total
shortwave can be converted from narrowband reflectance
using the formula developed earlier [Liang, 2001].
6. Summary
[52] High resolution incident PAR over land is needed by
ecological and hydrological modeling, but has not been
routinely generated from satellite observations. A new
method for estimating incident PAR from MODIS imagery
using the look-up table method has been presented above.
The unique feature of this new method is that both surface
reflectance and atmospheric properties are estimated simul-
taneously, while most existing algorithms require ancillary
information from other sources about aerosols and clouds.
The use of the look-up table also can enable us to avoid use
of the two-stream approximations in calculating multiple
scattering as used in most current algorithms since they may
result in significant errors. The table search is also very
expedient, so this method can be easily implemented for
regional and global applications. Its limitation is that it
cannot include many input variables. The linear interpola-
tion in the table search may introduce errors particularly
when the intervals in the table are large.
[53] The surface reflectances are estimated from the
temporal signatures using the minimum blue reflectance.
The look-up tables are established by using the MODTRAN
radiative transfer package. The method can separate direct
and diffuse components of incident PAR, and the total
?9?
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
11 of 13
D15208
shortwave radiative flux (insolation) can be also determined
from the same look-up table.
[54] The method has been validated using FLUXNET
observations from seven stations. The validation results
indicate that this method is reasonably accurate. Although
more validation activities are on the way, we feel confident
at this point that this method is very reliable for mapping
incident PAR from MODIS data.
[55] Because of the importance of diffuse PAR in
modeling and other applications, validating the diffuse
PAR product is critical. However, such ground measure-
ments are rare. Greater efforts should be made to measure
diffuse PAR component in the current or future observa-
tion networks.
[56] MODIS estimates of incident PAR are instantaneous;
the ideas for calculating daily averages are also discussed.
The current algorithms for generating PAR products have
not usually considered topography. This issue is particularly
relevant for carbon cycle modeling when the spatial reso-
lution increases, since most forests are distributed over the
mountainous regions. These issues will be addressed in the
future. We plan to evaluate the impacts of the improved
incident PAR product on the MODIS GPP/NPP products,
and distribute this product through the University of Mary-
land Global Land Lover Facilities.
[57] Acknowledgments. The authors would like to thank the
FLUXNET program and the investigators of seven stations where the in
situ measurements were providedfor our validation, John Townshendfor his
comments on the early draft, and anonymous reviewers for their valuable
comments that have greatly improved the presentation of this paper. This
study is partially funded by NASA grants NNG05GD10G and NAG512892.
R. Liu was supported in part by the National Natural Science Foundation
from China (40471098).
References
Alados, I., I. Foyo-Moreno, and L. Alados-Arboledas (1996), Photosyntheti-
cally active radiation: measurements and modelling, Agric. For. Meteorol.,
78,121?131.
Berk, A., L. S. Bernstein, G. P. Anderson, P. K. Acharya, D. C. Robertson,
J. H. Chetwynd, and S. M. Adler-Golden (1998), MODTRAN cloud and
multiple scattering upgrades with application to AVIRIS, Remote Sens.
Environ., 65, 367?375.
Bonan, G. B., K. W. Oleson, M. Vertenstein, S. Levis, X. Zeng, Y. Dai, R. E.
Dickinson, and Z. L. Yang (2002), The land surface climatology of the
community land model coupled to the NCAR community climate model,
J. Clim., 15, 3123?3149.
Carder, K., R. Chen, and S. Hawes (2003), Instantaneous photosyntheti-
cally available radiation and absorbed radiation by phytoplankton,
MODIS Ocean Science Team Algorithm Theoretical Basis Document,
ATBD MOD-20, version 7, 24 pp., NASA Goddard Space Flight Cent.,
Greenbelt, Md., Feb. (Available at http://modis.gsfc.nasa.gov/data/atbd/
ocean_atbd.php).
Chandrasekhar, S. (1960), Radiative Transfer, Dover, Mineola, N. Y.
Dai, Y., et al. (2003),The CommonLand Model (CLM), Bull. Am.Meteorol.
Soc., 84(8), 1013?1023.
Dickinson, R. E. (1995), Land-atmosphere interaction, Rev. Geophys,
33(S1), 917?922.
Dye, D. (2004), Spectral composition and quanta-to-energy ratio of dif-
fuse photosynthetically active radiation under diverse cloud conditions,
J. Geophys. Res., 109, D10203, doi:10.1029/2003JD004251.
Dye, D. G., and R. Shibasaki (1995), Intercomparison of global PAR data
sets, Geophys. Res. Lett., 22, 2013?2016.
Eck, T. F., and D. G. Dye (1991), Satellite estimation of incident photo-
synthetically active radiation using ultraviolet reflectance, Remote Sens.
Environ., 38, 135?146.
Field, C. B., J. T. Randerson, and C. M. Malmstrom (1995), Global net
primary production: Combining ecology and remote sensing, Remote
Sens. Environ., 51, 74?88.
Foley, J. A., I. C. Prentice, N. Ramunkutty, S. Levis, D. Pollard, S. Sitch,
and A. Haxeltine (1996), An integrated biosphere model of land surface
Figure 12. Predicting daily average PAR from one or two
simultaneous PAR at MODIS overpass times using (a) one
afternoon MODIS PAR, (b) one morning MODIS PAR, and
(c) both.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
12 of 13
D15208
processes, terrestrial biosphere model of ecosystem dynamics, Global
Biogeochem. Cycles, 10, 603?628.
Frouin, R., and R. Pinker (1995), Estimating photosynthetically active ra-
diation (PAR) at the Earth?s surface from satellite observations, Remote
Sens. Environ., 51, 98?107.
Frouin, R., B. Franz, and M. Wang (2000), Algorithm to estimate PAR from
SeaWiFS data, NASA SeaWiFs technical report, version 1.2, 15 pp.,
NASA Goddard Space Flight Cent., Greenbelt, Md. (Available at
http://daac.gsfc.nasa.gov/oceancolor/documentation/OB_Documenta-
tion.shtml).
Gu, L., D. Baldocchi, S. B. Verma, T. A. Black, T. Vesala, E. M. Falge, and
P. R. Dowty (2002), Advantages of diffuse radiation for terrestrial eco-
system productivity, J. Geophys. Res., 107(D6), 4050, doi:10.1029/
2001JD001242.
Gu, L. H., D. D. Baldocchi, S. C. Wofsy, J. W. Munger, J. J. Michalsky, S. P.
Urbanski, and T. A. Boden (2003), Response of a deciduous forest to the
Mount Pinatubo eruption: Enhanced photosynthesis, Science, 299(5615),
2035?2038.
Gu, J., et al. (2004), Modeling carbon sequestration over the large-scale
Amazon basin, aided by satellite observations. Part I: Wet- and dry-
season surface radiation budget flux and precipitation variability based
on GOES retrievals, J. Appl. Meteorol., 43, 870?886.
Houghton, R. A., E. A. Davidson, and G. M. Woodwell (1998), Missing
sinks, feedbacks, and understanding the role of terrestrial ecosystems in
the global carbon balance, Global Biogeochem. Cycles, 12(1), 25?34.
Jacovides, C. P., F. S. Tymvios, D. N. Asimakopoulos, K. M. Theofilou,
and S. Pashiardes (2003), Global photosynthetically active radiation and
its relationship with global solar radiation in the eastern Mediterranean
basin, Theor. Appl. Climatol., 74(3?4), 227?233.
Kaufman, Y. J., D. Tanre?, L. Remer, E. F. Vermote, A. Chu, and B. N.
Holben (1997a), Operational remote sensing of tropospheric aerosol over
the land from EOS-MODIS., J. Geophys. Res., 102(D14), 17,051?
17,068.
Kaufman, Y. J., A. Wald, L. A. Lorraine, B. C. Gao, R. R. Li, and L. Flynn
(1997b), Remote sensing of aerosol over the continents with the aid of a
2.2 um channel, IEEE Trans. Geosci. Remote Sens., 35, 1286?1298.
Liang, S. (2001), Narrowband to broadband conversions of land surface
albedo, Remote Sens. Environ., 76, 213?238.
Liang, S. (2004), Quantitative Remote Sensing of Land Surfaces, 534 pp.,
John Wiley, Hoboken, N. J.
Monteith, J. L. (1972), Solar radiation and productivity in tropical eco-
systems, J. Appl. Ecol., 9, 747?766.
Monteith, J. L. (1977), Climate and the efficiency of crop production in
Britain, Philos. Trans. R. Soc. London, Ser. B, 281, 277?294.
Myneni, R. B., et al. (2002), Global products of vegetation leaf area and
fraction absorbed PAR from year one of MODIS data, Remote Sens.
Environ., 83(1?2), 214?231.
Pinker, R. T., and I. Laszlo (1992), Modeling surface solar irradiance for
satellite applications on a global scale, J. Appl. Meteorol., 31, 194?211.
Pinker, R. T., et al. (2003), Surface radiation budgets in support
of the GEWEX Continental-Scale International Project (GCIP) and the
GEWEX Americas Prediction Project (GAPP), including the North
American Land Data Assimilation System (NLDAS) project, J. Geophys.
Res., 108(D22), 8844, doi:10.1029/2002JD003301.
Potter, C. S., et al. (1993), Terrestrial ecosystem production: a process
model based on global satellite and surface data, Global Biogeochem.
Cycles, 7, 811?841.
Prentice, I. C., W. Cramer, S. P. Harrison, R. Leemans, R. A. Monserud, and
A. M. Solomon (1992), A global biome model based on plant physiology
and dominance, soil properties and climate, J. Biogeogr., 19, 117?134.
Prince, S. D. (1991), A model of regional primary production for use with
coarse resolution satellite data, Int. J. Remote Sens., 12, 1313?1330.
Prince, S. D., and S. N. Goward (1995), Global primary production: a
remote sensing approach, J. Biogeogr., 22(4?5), 815?835.
Rousseeuw, P. J. (1984), Least median of squares regression, J. Am. Stat.
Assoc., 79, 871?881.
Ruimy, A., L. Kergoat, and A. Bondeau (1999), Comparing global models
of terrestrial net primary productivity (NPP): Analysis of differences in
light absorption and light-use efficiency, Global Change Biol., 5, 56?64.
Running, S., and S. T. Gower (1991), Forest-BGC, a general model of
forest ecosystem processes for regional applications, II. Dynamic alloca-
tion and nitrogen budgets, Tree Physiol., 9, 147?160.
Running, S. W., R. Nemani, J. M. Glassy, and P. Thornton (1999), MODIS
daily photosynthesis (PSN) and annual net primary production (NPP)
product (MOD17), MODIS Land Science Team Algorithm Theoretical
Basis Document (ATBD), ATBD-MOD16, version 3.0, 59 pp., NASA
Goddard Space Flight Cent., Greenbelt, Md., 29 Apr. (Available at http://
modis.gsfc.nasa.gov/data/atbd/land_atbd.php).
Running, S. W., P. E. Thornton, R. Nemani, and J. M. Glassy (2000),
Global terrestrial gross and net primary productivity from the Earth Ob-
serving System, in Methods in Ecosystem Science, edited by O. E. Sala et
al., pp. 44?57, Springer, New York.
Running, S. W., R. R. Nemani, F. A. Heinsch, M. S. Zhao, M. Reeves, and
H. Hashimoto (2004), A continuous satellite-derived measure of global
terrestrial primary production, Bioscience, 54(6), 547?560.
Salomonson, V., W. Barnes, P. Maymon, H. Montgomery, and H. Ostrow
(1989), MODIS: Advanced facility instrument for studies of the Earth as
a system, IEEE Trans. Geosci. Remote Sens., 27, 145?153.
Sellers, P. J., S. O. Los, C. J. Tucker, C. O. Justice, D. A. Dazlich, G. J.
Collatz, and D. A. Randall (1996), A revised land surface parameteriza-
tion (SiB2) for Atmospheric GCMs. Part II: The generation of global
fields of terrestrial biophysical parameters from satellite data, J. Clim.,
9, 706?737.
Stamnes, K., S. C. Tsay, W. Wiscombe, and K. Jayaweera (1988), Numeri-
cally stable algorithm for discrete-ordinate-method radiative transfer in
multiple scattering and emitting layered media, Appl. Opt., 27, 2502?
2509.
Tarpley, J. D. (1979), Estimating incident solar radiation at the surface from
geostationary satellite data, J. Appl. Meteorol., 18, 1172?1181.
Van Laake, P. E., and G. A. Sanchez-Azofeifa (2004), Simplified atmo-
spheric radiative transfer modelling for estimating incident PAR using
MODIS atmosphere products, Remote Sens. Environ., 91, 98?113.
Van Laake, P. E., and G. A. Sanchez-Azofeifa (2005), Mapping PAR using
MODIS atmosphere products, Remote Sens. Environ., 94(4), 554?563.
C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0C0
H. Fang, S. Liang, and T. Zheng, Department of Geography, University
of Maryland, College Park, MD 20742, USA. (sliang@geog.umd.edu)
R. Liu, Institute for Geographical Science and Natural Resources Re-
search, Chinese Academy of Sciences, Beijing 100101, China.
S. Running, School of Forestry, University of Montana, Missoula, MT
59812, USA.
S.-C. Tsay, Climate and Radiation Branch, NASA/Goddard Space Flight
Center, Greenbelt, MD 20771, USA.
D15208 LIANG ET AL.: PAR ESTIMATION FROM MODIS
13 of 13
D15208