remote sensing

Article
Mapping Maize Water Stress Based on UAV
Multispectral Remote Sensing
Liyuan Zhang 1,2 , Huihui Zhang 3 , Yaxiao Niu 1,2 and Wenting Han 1,4, *
1 College of Mechanical and Electronic Engineering, Northwest A&F University, Yangling 712100, Shaanxi,
China; [email protected] (L.Z.); [email protected] (Y.N.)
2 Key Laboratory of Agricultural Internet of Things, Ministry of Agriculture, Yangling 712100, Shaanxi, China
3 Water Management and Systems Research Unit, USDA-ARS, 2150 Centre Avenue, Bldg. D.,
Fort Collins, CO 80526, USA; [email protected]
4 Institute of Soil and Water Conservation, Northwest A&F University, Yangling 712100, Shaanxi, China
* Correspondence: [email protected]; Tel.: +86-029-8709-1325




Received: 30 January 2019; Accepted: 9 March 2019; Published: 13 March 2019


Abstract: Mapping maize water stress status and monitoring its spatial variability at a farm scale
are a prerequisite for precision irrigation. High-resolution multispectral images acquired from an
unmanned aerial vehicle (UAV) were used to evaluate the applicability of the data in mapping water
stress status of maize under different levels of deficit irrigation at the late vegetative, reproductive
and maturation growth stages. Canopy temperature, field air temperature and relative humidity
obtained by a handheld infrared thermometer and a portable air temperature/relative humidity
meter were used to establish a crop water stress index (CWSI) empirical model under the weather
conditions in Ordos, Inner Mongolia, China. Nine vegetation indices (VIs) related to crop water stress
were derived from the UAV multispectral imagery and used to establish CWSI inversion models.
The results showed that non-water-stressed baseline had significant difference in the reproductive
and maturation stages with an increase of 2.1 ◦ C, however, the non-transpiring baseline did not
change significantly with an increase of 0.1 ◦ C. The ratio of transformed chlorophyll absorption in
reflectance index (TCARI) and renormalized difference vegetation index (RDVI), and the TCARI and
soil-adjusted vegetation index (SAVI) had the best correlations with CWSI. R2 values were 0.47 and
0.50 for TCARI/RDVI and TCARI/SAVI at the reproductive and maturation stages, respectively;
and 0.81 and 0.80 for TCARI/RDVI and TCARI/SAVI at the late reproductive and maturation stages,
respectively. Compared to CWSI calculated by on-site measurements, CWSI values retrieved by
VI-CWSI regression models established in this study had more abilities to assess the field variability
of crop and soil. This study demonstrates the potentiality of using high-resolution UAV multispectral
imagery to map maize water stress.

Keywords: crop water stress index (CWSI); vegetation index; regression model; non-water-stressed
baseline; non-transpiring baseline

1. Introduction
Water resource scarcity is one of the most important challenges for agricultural water management,
especially in arid and semi-arid areas all over the world. Due to drought associated with climate change,
agricultural water resources will be reduced, it is necessary to achieve maximum production per unit
of applied irrigation water. As one of the important ways to enhance water-use efficiency, regulated
deficit irrigation may solve this dilemma [1–4]. However, to achieve a delicate balance between yield
and irrigation water, effective monitoring methods for crop water stress are necessary [5–7].

Remote Sens. 2019, 11, 605; doi:10.3390/rs11060605 www.mdpi.com/journal/remotesensing

Remote Sens. 2019, 11, 605 2 of 24

Currently, crop water stress can be detected based on soil moisture content, crop physiological
characteristics (e.g., stomatal conductance, leaf water potential) and remote-sensing technology [8].
On-site measurements of soil water content and crop physiological characteristics are time-consuming,
laborious and costly, and do not take into account the spatial variability of soil and crops [9,10].
Measurements of canopy temperature or canopy reflectance based on remote-sensing technology have
the advantages of being easy, non-destructive and of low labor intensity [8,11].
Since numerous data from various remote sensing imagery platforms have become available,
studies of monitoring crop water stress have been done by satellite [12–14] and aerial [15–18]
remote-sensing technology. However, crop monitoring applications are highly demanding with respect
to temporal frequency (1–3 days revisit) of key image acquisition, and spatial resolution requirements
are driven by the minimum management unit that can be in the range of 1–10 m [19]. Satellite
remote sensing has the disadvantages of low spatial-temporal resolution and vulnerability to weather.
For instance, Landsat and SPOT has spatial resolution of 30 and 20 m, and revisit time 16 and 26 days,
respectively. It is often difficult to meet the monitoring need of crop water stress at a farm scale. Aerial
remote sensing has the disadvantages of high operational difficulty and high operating cost. Due to the
advantages of low cost, easy to construct, convenient transportation, high flexibility, short operating
cycle and high spatial-temporal resolution, an unmanned aerial vehicle (UAV) remote-sensing system
could collect crop information in a desirable spatial-temporal resolution [20–22]. This makes it more
suitable to monitor crop water stress quickly and meticulously at a farm scale.
Due to the principle that crop canopy temperature (Tc) rises with increasing water stress in
comparison to crop under well irrigation condition [23], several researchers integrated thermal
cameras into UAV platforms to evaluate crop water stress of the entire field. Zhang et al. [24] used
canopy temperature characteristics derived from UAV thermal images to monitor cotton water status.
The results suggested canopy temperature standard deviation (CTSD) was sensitive to crop water stress
with a coefficient of determination (R2 ) of 0.88 between CTSD and stomatal conductance. Martínez
et al. [25] assessed the potential of UAV thermal imaging for monitoring plant water stress in a
commercial sugar beet field. The results indicated that a lightweight UAV canopy temperature system
including an infrared thermometer and a microcontroller board as data logger, was robust and reliable,
but there were some constraints related to weather conditions and delimitation of the area covered by
the infrared sensor.
However, Tc is not only affected by water supply but also by the actual micro-meteorological
conditions (such as, air temperature (Ta), relative humidity (RH) and vapor pressure deficit (VPD)).
Therefore, aiming to reduce the influence of meteorological conditions, several temperature-based
indices were developed over the last few decades [17]. The crop water stress index (CWSI) based on
canopy temperature and some meteorological parameters has been investigated as an effective method
to monitor the water status of various plants by many previous researchers, such as maize [26–29],
cotton [30,31], grapevine [16,32–34], peaches [35–37] and olives [38–40]. Good relationships between
CWSI and other on-site crop water stress indicators have been shown, such as soil water content [41–44],
leaf water potential [16,45,46] and stomatal conductance [32,33,47]. There are two widely used CWSI
models: the empirical model proposed by Idso et al. [48] and the theoretical model proposed by
Jackson et al. [49]. Compared to the theoretical model, the empirical model is easier to establish and
only needs Tc, Ta and RH.
When establishing the CWSI empirical model, there are two common methods to establish
non-water-stressed baseline (NWSB) and non-transpiring (stomata fully closed) baseline (NTB), namely
using wet and dry reference, and using the linear segment in the plot of dTLL (the lower limit of
canopy-air temperature difference) vs. VPD (vapor pressure deficit). For example, Wang et al. [44]
established the CWSI empirical model for winter wheat planted in Luancheng area, Hebei Province,
China, by using wet and dry reference to establish NWSB and NTB. However, to ensure the monitoring
effect of crop water stress, it is necessary to select the wet and dry reference position in practical
applications carefully. Taghvaeian et al. [42] established the CWSI empirical model for maize planted
Remote Sens. 2019, 11, 605 3 of 24

in Northeastern Colorado, USA by using the linear segment in the plot of dTLL vs. VPD. Compared to
using wet and dry reference, it is easier to establish NWSB and NTB by using the linear segment in the
plot of dTLL vs. VPD in practical applications.
As the key parameter to estimate CWSI, Tc conventionally relies on manual or continuous
point measurements [4,31,50]. However, the ground-based handheld thermometer approach is
labor-intensive, costly, and impractical when monitoring large areas; furthermore, water stress status
at the agricultural fields commonly has spatial variability due to spatially heterogeneous biophysical
conditions between plants or even within a plant [51–54]. To monitor crop water stress and its
spatial variability at a field scale, several researchers started to obtain CWSI maps based on a UAV
remote-sensing system. Bellvert et al. [55] characterized the spatial variability in water status across
vineyards by CWSI maps. The results showed that the CWSI maps thus obtained by UAV thermal
imagery were effective in assessing the spatial variability of water stress across the vineyard. Zhang et
al. [30] evaluated cotton water stress and its spatial variability by UAV thermal remote-sensing system.
The results showed that the R2 between CWSI and leaf stomatal conductance of cotton planted in
Yangling, Shaanxi, China was 0.84.
However, due to the load capacity limitation of UAV platforms, the thermal infrared sensors
currently mounted are mostly lightweight uncooled thermal infrared cameras that have the problems
of unstable canopy temperature acquisition, low mosaic precision (e.g., number of tie-points, number
of geolocated images and mean re-projection error etc.) and needing complex pretreatments (e.g.,
radiometric calibration, temperature correction and canopy temperature extraction) [17,56]. Compared
to the thermal infrared system, a UAV-based multispectral remote-sensing system has the advantages
of stable information acquisition and mature mosaic technology. It could be a supplement to
ground-based infrared thermometry to obtain CWSI maps with a high spatial-temporal resolution at
the field scale.
A vegetation index (VI), which is the mathematical operation result of reflectivity of vegetation
and other ground objects in the corresponding band [57], can be used to monitor surface vegetation
condition. VI has been successfully applied to land cover classification, environmental change, and
drought monitoring [33]. Water stress could affect crop canopy structure, photosynthesis and leaf
moisture content. VIs used to monitor crop water stress can be divided into four categories: structure
index, chlorophyll index, lutein index and moisture index [8]. The structure and chlorophyll indices
based on visible, near infrared and red-edge bands are more widely used, such as normalized difference
vegetation index (NDVI) [58], renormalized difference vegetation index (RDVI) [33], soil-adjusted
vegetation index (SAVI) [59], optimization of soil-adjusted vegetation index (OSAVI) [60], transformed
chlorophyll absorption in reflectance index (TCARI) [60]. To reduce the sensitivity of TCARI to canopy
structure changes, TCARI/OSAVI was used by Haboudane et al. [60]. When structure index was used
to describe canopy structure, its effect would be affected by the soil, crops, and planting structure.
To reduce the sensitivity of TCARI to canopy structure changes, in this study NDVI, RDVI, SAVI and
OSAVI were used, respectively.
To map maize water status and monitor its spatial variability at a farm scale, in this study,
the CWSI empirical model was developed by the ground-based data and high-resolution UAV-based
multispectral images of maize under various levels of deficit irrigation at different growth stages.
The linear segment in the plot of dTLL vs. VPD was used to establish the NWSB and NTB. VIs related
to water stress derived from multispectral imagery were used to establish CWSI inversion model.
Since, the local conditions, irrigation technique, and management practices can greatly influence the
applicability of remote sensing systems for water status evaluation. The specific objectives of this
study were:

(1) To develop the NWSB and NTB for maize under the weather condition in Ordos, Inner
Mongolia, China;
(2) To evaluate UAV-based multispectral VIs that are sensitive to maize water stress and establish
CWSI regression models;
Remote Sens. 2019, 11, 605 4 of 24

(3) To obtain CWSI maps with high spatial-temporal resolution at the field scale using the CWSI
regression models.

2. Materials and Methods

2.1. Study Site

This study was conducted in a 1.13 ha research field (40◦ 260 0.29”N, 109◦ 360 25.99”E, Elev. 1010 m),
located in Zhaojun Town, Dalate Banner, Ordos, Inner Mongolia, China. Soil samples at the depths
of 30, 60 and 90 cm were collected from the sampling plot (yellow circle in Figure 1a) to estimate
the field hold capacity (volumetric) which were 13.6%, 13.2% and 12.3%, respectively. The average
permanent wilting point of the soil profile was assumed to be 5.6% (volumetric) according to Lv et
al. [61]. The soil bulk density was 1.56. The mixture of soil samples at the depths of 30, 60 and 90 cm
were used to estimate other soil characteristics. The soil type was a loamy sand (80.7% sand, 13.7%
powder, and 5.6% clay), using the United States Department of Agriculture soil taxonomy. The soil pH,
organic matter and C content were 9.27, 47.17 g/Kg and 27.35 g/Kg, respectively. Maize (Junkai 918)
was planted on 20 May 2017, with 0.58 m row spacing and 0.25 m plant spacing and the row direction
was East-West. The maize emerged on June 1, headed on July 20, was harvested on September 7 (silage)
with
RemoteaSens.
110-day lifespan.
2019, 11, x FOR PEER REVIEW 5 of 25

(a) (b)
Figure 1.
Figure 1. Aerial
Aerial view
view of
of the
the experimental
experimental field
field indicating
indicating treatment
treatment region
region division
division (a),
(a), the
the location
location
of sampling plots (a) and sampling sites (b). Air temperature and humidity, and canopy temperature
of sampling plots (a) and sampling sites (b). Air temperature and humidity, and canopy temperature
were taken
were taken at
at each
each sampling
sampling site.
site.

2.2. Experimental Design

Table 1. Experimental treatments and total applied water depth (percentage of full irrigation
treatment in parentheses) that includes amount of irrigation and precipitation in the late vegetative,
The study field was divided into five treatment (TRT) regions (Figure 1a) with different levels of
reproductive and maturation stages in 2017.
deficit irrigation at the crop late vegetative, reproductive and maturation growth stages. A 12 × 12 m2
area within each TRT region was chosen as a Appliedsampling plot Depth/mm
Water for the TRT (Figure 1a) and data collection
from five sampling Treatment
sites withinLate Vegetativeplot
the sampling Reproductive
were indicated Maturation
in Figure 1b.Total
The non-stressed crop
evapotranspiration (ET) was estimated(07.04–07.28) (07.29–08.20)
by reference (08.21–09.07)and single crop coefficient
evapotranspiration
approach [62]. The crop TRT 1coefficient
188 was
(100%)
equal to132 (100%)
0.55, 1.09 and820.45
(100%) 402
in the initial, mid-season and
TRT 2 158 (84%) 128 (97%) 43 (52%) 329
late-season development stages. The total crop water requirement of the full watered maize during
TRT 3 158 (84%) 125 (95%) 43 (52%) 326
the late vegetation, reproductive and maturation stages was 407 mm which was close to the total
TRT 4 158 (84%) 91 (69%) 23 (28%) 272
TRT 5 158 (84%) 124 (94%) 82 (100%) 365

2.3. Measurement of Maize Canopy Temperature

Maize canopy temperature (Tc) and UAV multispectral data were synchronously collected
under clear sky conditions between 11:00–13:00 local time (Chinese standard time, 11:44–13:44)
Remote Sens. 2019, 11, 605 5 of 24

applied water (402 mm) of the control treatment region (TRT 1). Therefore, different levels of deficit
irrigation were described by the percentages of applied water depth of TRT 1 during the late vegetation,
reproductive and maturation stages. For example, during the late vegetation stage, 84% of the applied
water depth at TRT 1 was applied to TRT 2 (Table 1)

Table 1. Experimental treatments and total applied water depth (percentage of full irrigation treatment
in parentheses) that includes amount of irrigation and precipitation in the late vegetative, reproductive
and maturation stages in 2017.

Applied Water Depth/mm

Treatment Late Vegetative Reproductive Maturation
Total
(07.04–07.28) (07.29–08.20) (08.21–09.07)
TRT 1 188 (100%) 132 (100%) 82 (100%) 402
TRT 2 158 (84%) 128 (97%) 43 (52%) 329
TRT 3 158 (84%) 125 (95%) 43 (52%) 326
TRT 4 158 (84%) 91 (69%) 23 (28%) 272
TRT 5 158 (84%) 124 (94%) 82 (100%) 365

Irrigation was applied during the growing season by using a center pivot sprinkler system
(Valmont Industries, Inc., Omaha, NE, USA) that has two spans and one end gun with a total length
of 143.7 m. More detailed information about the center pivot sprinkler system can be found in Li et
al. [63]. The uniformity test of water application of the center pivot irrigation system was conducted in
accordance with the standard ANSI/ASAE S436.1 and ISO 11545. The coefficient of uniformity of the
first span (research field) calculated by the modified formula of Heermann and Hein [64] was 82.7%
under a walking speed of 20% of the full speed and 88.3% under a walking speed of 40% of the full
speed using R3000 sprinklers. The amount of water applied to each treatment was measured and
recorded by MIK-2000H flow meters (Meacon Automation Technology Co., Ltd, Hangzhou, China).
The actual amount of irrigation and rainfall at each grown stage were shown in Table 1. In order
to eliminate the interference of nutritional stress and weeds, fertilizer and herbicide were applied
according to the planting experience.

2.3. Measurement of Maize Canopy Temperature

Maize canopy temperature (Tc) and UAV multispectral data were synchronously collected under
clear sky conditions between 11:00–13:00 local time (Chinese standard time, 11:44–13:44) during the
study period (2017.08.06~29, including the reproductive and maturation stages). Tc was measured by a
handheld infrared thermometer (RAYTEK, ST60+, Raytek Inc, Santa Cruz, CA, USA). The temperature
range and spectral range were −32–600 ◦ C and 8–14 µm, respectively. The measurement error was
±1% of reading or ±1 ◦ C, the lager value was adopted in the practical application. The emissivity
value was set to 0.97. To avoid the interference of the soil, the handheld infrared thermometer was
used to sweep the canopy (about 120◦ ) perpendicular to the row with 30 cm above the canopy and
15◦ horizontal angle (Figure 2). At each sampling site three measurements of Tc were taken and the
average value of these three readings was used to represent the sampling site. The average value of Tc
measured at the five sampling sites (yellow circle in Figure 1b) was used to represent the status of the
sampling plot. During each data collection day, a total of 15 canopy temperature measurements were
obtained for each sampling plot. During the reproductive and maturation stages, the average plant
heights were 263, 258, 252, 253 and 261 cm for TRTs 1~5, respectively.
Remote Sens. 2019, 11, 605 6 of 24
Remote Sens. 2019, 11, x FOR PEER REVIEW 6 of 25

Figure
Figure 2.
2. Schematic
Schematic indicating
indicating measurement method of
measurement method of canopy
canopy temperature
temperature in
in each
each sampling
sampling site.
site.

2.4. Meteorological Data

2.4. Meteorological Data
The field
fieldair
airtemperature
temperature andandrelative humidity
relative (Ta_f(Ta_f
humidity and RH_f) of eachof
and RH_f) sample
each plots wasplots
sample obtained
was
by a handheld
obtained hygrometer
by a handheld (Klimalogg
hygrometer Pro, TFA, Pro,
(Klimalogg Germany). During theDuring
TFA, Germany). study period Ta_f
the study and RH_f
period Ta_f
were measured
and RH_f synchronously
were measured with Tc atwith
synchronously a height
Tc atof 2.2 m above
a height of 2.2 mtheabove
ground. The average
the ground. of Ta_f
The average
and
of Ta_f and RH_f collected at the five sampling sites (Figure 1b) represented the status of plot.
RH_f collected at the five sampling sites (Figure 1b) represented the status of one sampling one
The weather
sampling data
plot. Thewere measured
weather at anmeasured
data were automated at weather stationweather
an automated adjacentstation
to theadjacent
researchtofield,
the
which included rainfall, Ta and RH, net solar radiation, and wind speed (2 m
research field, which included rainfall, Ta and RH, net solar radiation, and wind speed (2 m above above the reference
grass surface).
the reference In addition
grass to addition
surface). In rainfall, to
therainfall,
data acquisition interval was
the data acquisition 30 min.
interval was 30The main
min. Themean
main
meteorological data during
mean meteorological data the study
during period
the studyincluding reproductive
period including stage (07.29–08.20)
reproductive and maturation
stage (07.29–08.20) and
stage (08.21–29) were shown in Table 2. The changing curves of meteorological
maturation stage (08.21–29) were shown in Table 2. The changing curves of meteorological data were data were shown in
Figure
shown A1.
in Figure A1.

Table 2. The
Table 2. main mean
The main mean meteorological
meteorological data during the
data during the study
study period
period including
including the
the reproductive
reproductive stage
stage
and maturation stages in 2017.
and maturation stages in 2017.

Reproductive
Reproductive Maturation
Maturation
Parameter
Parameter (07.29–08.20)
(07.29–08.20) (08.21–29)
(08.21–29)
Mean Mean
air temp./ ◦
air temp./°C
C 22.1122.11 17.21
17.21
Max.Max.
air temp./ ◦C
air temp./°C 31.3131.31 25.46
25.46
Min.Min.
air temp./ ◦C 13.6113.61 9.24
air temp./°C 9.24
Min.Min.
relative humidity/%
relative humidity/% 29.7829.78 33.23
33.23
Mean
Mean net solarnet solar
radiation/MJ·m −2·day−1
−2 ·day−1 10.9810.98 3.00
3.00
radiation/MJ · m
Mean wind speed/m·s −1 0.47 0.28
Mean wind speed/m·s−1 0.47 0.28

2.5. Soil Water Content Data

2.5. Soil Water Content Data
Soil water content (SWC) was measured two or three times each week on the days before or after
irrigationwater
Soil in each content (SWC)
sampling plotwas measured
by the twogravimetric
traditional or three times each At
method. weekeachonsampling
the daysplot,
before or
three
after irrigation in each sampling plot by the traditional gravimetric method. At each
sampling locations were randomly chosen around the center. At each sampling site, soil samples at sampling plot,
three sampling
the depths of 30,locations
60 and 90were randomly
cm were chosen
collected by around theaugers.
using soil center.Soil
At each sampling
samples site,in
were put soil samples
aluminum
at the depths of 30, 60 and 90 cm were collected by using soil augers. Soil samples
boxes to avoid the influence of evaporation. An electronic balance (DJ-2000A, Suzhou Jiangdong were put in
aluminum boxes to avoid the influence of evaporation. An electronic balance (DJ-2000A,
Precision Instrument CO., Ltd., Suzhou, China) with a measurement accuracy of 0.01 g was used to Suzhou
Jiangdong
measure the Precision
weightInstrument CO., Ltd., boxes
of the aluminum Suzhou, China)
with andwith a measurement
without soil samples.accuracy
Afterof weighting,
0.01 g was
used to measure the weight of the aluminum boxes with and without soil samples.
uncapped aluminum boxes were dried to a constant weight by putting in an oven under 105 °C. SWC After weighting,
(volumetric) could be calculated by using Equation (1):

W1 − W2
SWC = ×1.56 (1)
W2 − W3
Remote Sens. 2019, 11, 605 7 of 24

uncapped aluminum boxes were dried to a constant weight by putting in an oven under 105 ◦ C. SWC
(volumetric) could be calculated by using Equation (1):

W1 − W2
Remote Sens. 2019, 11, x FOR PEER REVIEW SWC = × 1.56 7 of(1)
25
W2 − W3

where W
where W11presents
presentsthe
theweight
weight of aluminum
of aluminum boxboxand and soil sample
soil sample beforebefore
drying;drying; W2 presents
W 2 presents the
the weight
weight
of of aluminum
aluminum box sample
box and soil and soil sample
after after
drying; W 3drying;
presents Wthe
3 presents thealuminum
weight of weight ofbox
aluminum
without box
soil
without1.56
sample; soil is
sample; 1.56 issoil
the specific thebulk
specific soil of
density bulk
thisdensity of this
study site. Thestudy site. The
changing changing
curve curve of
of the average the
SWC
average SWCat(volumetric)
(volumetric) the depths ofat30,the60depths
and 90ofcm30,during
60 andthe
90reproductive
cm during the andreproductive
maturation and maturation
stages is shown
stages
in is shown
Figure 3. in Figure 3.

Figure 3. The changing curve of the average soil water content (SWC, volumetric) at the depths of 30,
Figure 3. The changing curve of the average soil water content (SWC, volumetric) at the depths of 30,
60 and 90 cm during the reproductive and maturation stages in 2017. The black, red, green, blue and
60 and 90 cm during the reproductive and maturation stages in 2017. The black, red, green, blue and
purple solid lines represent treatment regions (TRTs) 1~5, respectively. The red dotted line is boundary
purple solid lines represent treatment regions (TRTs) 1~5, respectively. The red dotted line is
between the late vegetation and reproductive stages. The black dotted line is boundary between the
boundary between the late vegetation and reproductive stages. The black dotted line is boundary
reproductive and maturation stages. The blue dotted lines indicate dates when data was taken for the
between theofreproductive
calculation and maturation
non-water-stressed stages. The baselines
and non-transpiring blue dotted lines
of the indicate
crop waterdates
stresswhen
indexdata was
(CWSI).
taken for the calculation of non-water-stressed and non-transpiring baselines of the crop water stress
2.6. Unmanned Aerial Vehicle (UAV) Multispectral System and Data Collection
index (CWSI).

Three types of UAV remote-sensing systems are generally used to collect agricultural information,
2.6. Unmanned Aerial Vehicle (UAV) Multispectral System and Data Collection
namely fixed-wing, single-rotor and multi-rotors UAVs [65]. When selecting a kind of UAV platform,
Three
factors suchtypes of UAV
as takeoff remote-sensing
requirements, coveragesystems are capacity,
area, load generallywind usedresistance,
to collectandagricultural
cost and
information, namely
maintenance cost arefixed-wing,
needed to be single-rotor
considered and multi-rotors
[56]. In this study,UAVs [65]. WhenUAV
a hexacopter selecting a kind of
multispectral
UAV platform, factors
remote-sensing systemsuch as takeoff
(Figure 4) wasrequirements,
developed with coverage
a Pixhawkarea, autopilot
load capacity, wind
(CUAV, resistance,
Guangzhou,
and cost and maintenance cost are needed to be considered [56]. In this
China), a RedEdge multispectral camera (MicaSense, Inc., WA, USA) and a MOY brushless gimbal study, a hexacopter UAV
multispectral remote-sensing
(Moyouzhijia, Huizhou, China). system
This(Figure
system 4)
hadwas thedeveloped
advantages with a Pixhawk
of easy autopilot
to takeoff (CUAV,
and land, low
Guangzhou, China), a RedEdge multispectral camera (MicaSense, Inc., WA,
operation cost and maintenance cost. The core of the UAV system was a Flight Control board, Pixhawk, USA) and a MOY
brushless
based gimbal
on one (Moyouzhijia,
32 bit STM32F427Huizhou,
Cortex M4 China).
core This
with system had the
float point unitadvantages
and one 32ofbiteasy to takeoff
STM32F103
and land,
failsafe low operation
co-processor, cost and
which maintenance with
communicated cost. The
the core of the UAVcontrollers
six brushless system wasby a Flight Control
six two-wire
board, Pixhawk, based on one 32 bit STM32F427 Cortex M4 core with float point
bi-directional wires (pulse width modulation). Pixhawk integrated two gyroscopes, two accelerometers, unit and one 32 bit
STM32F103
one failsafe co-processor,
magnetometer, which communicated
and one barometer. with the six
The gimbal maintained thebrushless controllers
multispectral cameraby under
six two-a
wire bi-directional
relatively stable state wires
during(pulse width
the flight modulation).
to obtain good quality Pixhawk integrated
multispectral twoThe
images. gyroscopes, two
flight planning
accelerometers,
was conducted with one magnetometer,
ground control and onesoftware,
station barometer. The gimbal
Mission Planner, maintained the multispectral
which allowed the user to
camera under
generate a routea relatively stable
of waypoints as state duringofthe
a function theflight
sensorto field
obtain of good
view quality multispectral
(FOV), the images.
degree of overlap
between images and the ground resolution needed. Mission Planner could also display real-timewhich
The flight planning was conducted with ground control station software, Mission Planner, flight
allowed
data. the user
Its main to generate
technical a routewere
parameters of waypoints
shown inas a function
Table 3. of the sensor field of view (FOV), the
degree of overlap between images and the ground resolution needed. Mission Planner could also
display real-time flight data. Its main technical parameters were shown in Table 3.
Remote Sens. 2019, 11, 605 8 of 24

Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 25

Figure
Figure 4.
4. Schematic
Schematic indicating
indicatingthe
the main
main components
componentsof
of UAV
UAV multispectral
multispectral remote
remote sensing
sensing system
system
developed
developed in
in this
this study.
study.

Table 3. Main parameters of UAV multispectral image acquisition system.

Table 3. Main parameters of UAV multispectral image acquisition system.
Parameter
Parameter Value
Value
Wheelbase
Wheelbase 900 mm900 mm
TakeoffTakeoff
weight weight 6 kg 6 kg
Payload Payload 2 kg 2 kg
UAVUAV Flight time 18 min
Flight time 18 min
Communication
Communication radius radius 3 km 3 km
Speed 5 m/s
Speed 5 m/s
Camera MicaSense RedEdge
Camera MicaSense RedEdge
Imager resolution 1280 × 960 pixels
ImagerImager
resolution
size 4.8 mm1280
× ×3.6
960
mmpixels
Imager size 4.8 mm × 3.6 mm
RedEdge Blue (475 nm), Green (560 nm), Red (668 nm),
Spectral bands Blue (475 nm), Green (560 nm), Red (668
Near infrared (IR) (840 nm), Red-edge (717 nm)
SpectralLens
bands
focal length 5.5 mmnm), Near
fixed infrared (IR) (840 nm), Red-edge
lens
RedEdge
Lens field of view ◦ (717 nm)
47.2 horizontal field of view (HFOV)
Lens focal length
Weight 150 g5.5 mm fixed lens
Dimension
Lens field of view 121 mm × horizontal
47.2° 66 mm × 46 mmof view (HFOV)
field
Weight 150 g
Dimension 121 mm × 66 mm × 46 mm
During the study period (2017.06.20~08.29), eleven UAV flights were conducted on sunny day
between 11:00–13:00 local time (Chinese standard time, 11:44–13:44) with the RedEdge camera lens
During the study period (2017.06.20~08.29), eleven UAV flights were conducted on sunny day
vertically downwards, and 80% heading and side overlap. The flight height, speed and ground sample
between 11:00–13:00 local time (Chinese standard time, 11:44–13:44) with the RedEdge camera lens
distance were 70 m (relative flying height), 5 m/s and 4.7 cm, respectively; 2185 images (5 bands)
vertically downwards, and 80% heading and side overlap. The flight height, speed and ground
were collected during a single flight and Pix4DMapper software was used for image mosaicking.
sample distance were 70 m (relative flying height), 5 m/s and 4.7 cm, respectively; 2185 images (5
The multispectral images of diffuse reflector (reflectivity 58%, size 3 × 3 m, Group VIII, USA) collected
bands) were collected during a single flight and Pix4DMapper software was used for image
simultaneously at the same height were used to perform radiometric correction. Then, spectral
mosaicking. The multispectral images of diffuse reflector (reflectivity 58%, size 3 × 3 m, Group VIII,
reflectance of the objects was obtained.
USA) collected simultaneously at the same height were used to perform radiometric correction. Then,
spectral reflectance of the objects was obtained.

2.7. Vegetation Indices’ Selection and Calculation

To establish a regression mode between UAV-based multispectral VI and CWSI, nine VIs, NDVI,
RDVI, SAVI, OSAVI, TCARI, TCARI/NDVI, TCARI/RDVI, TCARI/SAVI, and TCARI/OSAVI were
selected. Their calculation formulas are as follows:
Remote Sens. 2019, 11, 605 9 of 24

2.7. Vegetation Indices’ Selection and Calculation

To establish a regression mode between UAV-based multispectral VI and CWSI, nine VIs, NDVI,
RDVI, SAVI, OSAVI, TCARI, TCARI/NDVI, TCARI/RDVI, TCARI/SAVI, and TCARI/OSAVI were
selected. Their calculation formulas are as follows:
ρnir − ρred
NDV I = (2)
ρnir + ρred

ρ − ρred
RDV I = √ nir (3)
ρnir + ρred
(1 + 0.5) ∗ (ρnir − ρred )
SAV I = (4)
ρnir − ρred + 0.5
(1 + 0.16) ∗ (ρnir − ρred )
OSAV I = (5)
ρnir − ρred + 0.16
h i
TCARI = 3 (ρrededge − ρred ) − 0.2(ρrededge − ρ green ) ∗ (ρrededge /ρred ) (6)

where ρnir , ρred , ρrededge and ρ green are reflectance values of ground objects in near-infrared, red,
red-edge and green bands. A linear regression analysis between CWSI and VI was conducted to
obtain the VI vs. CWSI regression model based on CWSI and VI of sampling plots. For statistical
analysis, R programming language (R-3.4.3, https://www.r-project.org/) and the lm() function
were used. The coefficient of determination (R2 ) and root mean square error (RMSE) were used
as evaluating indicators.

2.8. Crop Water Stress Index (CWSI) Calculation

The empirical CWSI could be calculated by using Equations (7)–(10):

dTm − dTLL
CWSI = (7)
dTUL − dTLL

where dTm , dTLL , dTUL are actual measurement, lower limit and upper limit of canopy-temperature
difference, respectively. The upper and lower limit canopy-air temperature difference are the so-called
non-water-stressed baseline (NWSB) and non-transpiring baseline (NTB). They could be calculated by
Equations (8) and (9) [22,37]:
dTLL = m × VPD + b (8)

dTUL = m × VPG + b (9)

where m, b are slope and intercept, respectively; VPD (vapor pressure deficit) is related to Ta and RH,
and calculated by the following Equation (10) with the description of Walter et al. [66]:

17.27Ta 100 − RH
VPD = 0.6108 × exp( )×( ) (10)
Ta + 237.3 100

where Ta is air temperature, RH is relative humidity. Vapor pressure gradient (VPG) is the change of
VPD when the Ta increase by b, and calculated by the following Equation (11):

17.27( Ta + b) 100 − RH
VPG = 0.6108 × exp( )×( ) − VPD (11)
( Ta + b) + 237.3 100

where b is the intercept in Equations (8) and (9). When CWSI is 0, it indicates no water stress,
while CWSI is 1, it indicates most severe stress.
After effective irrigation/rainfall, it is assumed that the soil water deficit is fully replenished and
crop could obtain sufficient water. Thus, non-water-stressed condition exist and Tc-Ta values could be
Remote Sens. 2019, 11, 605 10 of 24

regarded as the values of the lower limit of canopy-air temperature difference [29,42]. Graph of hourly
dT LL vs.Sens.
Remote VPD has11,axlinear
2019, segment
FOR PEER REVIEWthat extended for 2 h before and 2 h after the solar noon. This linear 10 of 25
segment could be extracted to obtain the coefficients of Equations (8) and (9) by using simple linear
simple linear
regression regression
[42,67]. [42,67].
In this study, In this
after study,
effective after effective irrigation/rainfall,
irrigation/rainfall, clear were
clear sky sunny days sky sunny
chosendays
to
were chosencollect
continuously to continuously
(local timecollect (local time
10:00~14:00, with10:00~14:00,
an interval ofwith an interval
15 min) ofand
Tc, Ta_f 15 min)
RH_fTc,in Ta_f and
the full
RH_f in region
irrigation the full(TRT1)
irrigation region
during (TRT1) during
the reproductive thematuration
and reproductive andDuring
stages. maturation stages. During
the reproductive the
stage,
9 reproductive
August and 16 stage,
August 9 August and 16
2017 (clear skyAugust 2017 day
sunny day, (clear
of sky
yearsunny
(DOY)day,221 day
and of year
228) were(DOY)
chosen221and
and
228) were chosen and 28 groups of samples (Tc, Ta_f and RH_f) were obtained. During
28 groups of samples (Tc, Ta_f and RH_f) were obtained. During the maturation stage, 24 August 2017 the maturation
stage,sky
(clear 24sunny
August, day,2017
DOY (clear
236)sky
wassunny
chosen day,
andDOY 236) was
13 groups chosen were
of samples and 13 groups of
obtained. samples
When were
drawing
obtained.
the dTLL ~VPD Whenscatterplot, the dTLL ~VPD
drawing three-step moving scatterplot,
average wasthree-step
used tomoving average was
reduce volatility used totoreduce
according Idso
(1981) [48],according
volatility making the fitting
to Idso curve
(1981) could
[48], represent
making the lower
the fitting curvelimit
couldofrepresent
Tc-Ta better underlimit
the lower different
of Tc-
saturated
Ta betterVPD.
under different saturated VPD.

3.3.Results
Results

3.1.
3.1.Non-Water-Stressed
Non-Water-StressedBaselines
Baselines(NWSBs)
(NWSBs)and andNon-Transpiring
Non-Transpiring Baselines
Baselines(NTBs)
(NTBs)of of
CWSI
CWSI
ItItwas
wasobserved
observedthatthatthe
theslope
slopeand
andintercept
interceptofofNWSB
NWSBchanged
changedininthethematuration
maturationstage.
stage.AAsimilar
similar
observation
observation was made by previous researchers among crops, such as sunflower [43,68],wheat
was made by previous researchers among crops, such as sunflower [43,68], wheat[69],
[69],
barley [69] and maize [70]. Hence, different NWSBs were developed in this study for
barley [69] and maize [70]. Hence, different NWSBs were developed in this study for the reproductive the reproductive
and
andmaturation
maturation stages.
stages.The dTLL
The dTvs. VPD scatterplots were shown in Figures 5a and 5b, respectively.
LL vs. VPD scatterplots were shown in Figure 5a and Figure 5b,
The corresponding NWSBs were
respectively. The corresponding NWSBs were shown as Equations
shown as(12) and (13).
Equations (12) The average
and (13). The NWSBs of the
average NWSBs
reproductive and maturation stages were − 11.0 ◦ C and −8.9 ◦ C, rose by 2.1 ◦ C. In calculating CWSI,
of the reproductive and maturation stages were −11.0 °C and −8.9 °C, rose by 2.1 °C. In calculating
the NTBthe
CWSI, is estimated using theusing
NTB is estimated samethecoefficients developed
same coefficients for NWSB,
developed forand usingand
NWSB, VPG instead
using VPGofinstead
VPD.
The change in NWSB coefficients from the reproductive to maturation stage resulted
of VPD. The change in NWSB coefficients from the reproductive to maturation stage resulted a small a small change in
NTB estimations. The averageThe NTBs of the NTBs
reproductive and maturation ◦ ◦
change in NTB estimations. average of the reproductive andstages were 0.7
maturation stagesC and
were0.80.7C,°C
with a rise of 0.1 ◦ C. Taghvaeian et al. [42,67] also showed a negligible variation when to estimate
and 0.8 °C, with a rise of 0.1 °C. Taghvaeian et al. [42,67] also showed a negligible variation when to
maize
estimateNTBmaize
in northern
NTB inColorado.
northern Colorado.
=LL−=2.64
dTLLdT −2.64
× ×VPD + 0.42
VPD + R22 ==0.89)
0.42 ((R 0.89) (12)
(12)
2
=LL−=3.35
dTLLdT −3.35
× ×VPD + 2.96
VPD + R2 ==0.98)
2.96 ((R 0.98) (13)
(13)

(a) (b)
Figure5. Maize
Figure 5. Maize non-water-stressed
non-water-stressed baselines
baselines (NWSBs)(NWSBs) developed
developed in during
in this study this study during the
the reproductive
(a)reproductive (a) and
and maturation maturation
stages stages (b) in 2017.
(b) in 2017.

3.2. CWSI of Maize

Finally, the CWSI empirical model was obtained by using Tc, Ta_f and RH_f collected in the
sampling plot. Figure 6 depicts the daily changes of CWSI and irrigation/rainfall events for the TRT
1 (Figure 6a) and TRT 4 (Figure 6b) during 2017.08.06~29. CWSI decreased after irrigation/rainfall,
reached a minimum around the third day and then slowly increased, responding well to
irrigation/rainfall events. When irrigation/rainfall was less, no more than 4 mm, CWSI did not change
Remote Sens. 2019, 11, 605 11 of 24

3.2. CWSI of Maize

Finally, the CWSI empirical model was obtained by using Tc, Ta_f and RH_f collected in
the sampling plot. Figure 6 depicts the daily changes of CWSI and irrigation/rainfall events
for the TRT 1 (Figure 6a) and TRT 4 (Figure 6b) during 2017.08.06~29. CWSI decreased after
irrigation/rainfall, reached a minimum around the third day and then slowly increased, responding
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 25
well to irrigation/rainfall events. When irrigation/rainfall was less, no more than 4 mm, CWSI did not
change significantly.
significantly. The average
The average CWSI forCWSIthefor the different
different levels levels of deficit
of deficit irrigation
irrigation in theinreproductive
the reproductive
and
and maturation stages had a clear numerical gradient (Table 4). For example,
maturation stages had a clear numerical gradient (Table 4). For example, the average CWSI the average CWSI was
was 0.11
0.11 0.27
and and 0.27 for TRT
for TRT 1 and
1 and TRTTRT
4 in4 the
in the reproductive
reproductive stage,
stage, respectively.The
respectively. Theaverage
averageCWSI
CWSIwas
was 0.03,
0.03,
0.14 and
0.14 and 0.32
0.32 for
for TRT
TRT 1,1, TRT
TRT 2,
2, and
and TRT
TRT 44 in
in the
the maturation
maturation stage,
stage, respectively.
respectively.

(a) (b)
Figure 6. CWSI changing
changing curves
curvesfor
forTRT
TRT11(a)(a)and
andTRT
TRT4 4(b)
(b)during
during 2017.08.06~29.
2017.08.06~29. The
The depths
depths (mm)
(mm) of
individual
of irrigation
individual (I) (I)
irrigation and precipitation
and (P)(P)
precipitation events were
events wereplotted as vertical
plotted bars.
as vertical bars.

Table 4. Mean CWSI for each deficit irrigation treatment during the reproductive and maturation
Table 4. Mean CWSI for each deficit irrigation treatment during the reproductive and maturation
stages in 2017.
stages in 2017.
Treatment
Treatment Reproductive (08.06–08.20)
Reproductive (08.06–08.20) Maturation
Maturation (08.21–08.29)
(08.21–08.29)
TRT TRT
1 1 0.07 0.03 0.03
TRT TRT
2 2 0.12 0.14 0.14
TRT TRT
3 3 0.09
0.09 0.12 0.12
TRT TRT
4 4 0.27
0.27 0.32 0.32
TRT 5 0.05 0.03
TRT 5 0.05 0.03

3.3. Vegetation
3.3. Vegetation Indices
Indices
Figure 77 shows
Figure shows thethe changes
changes of of NDVI,
NDVI, RDVI,
RDVI, SAVI,
SAVI, OSAVI,
OSAVI,TCARI,
TCARI,TCARI/NDVI,
TCARI/NDVI, TCARI/RDVI,
TCARI/RDVI,
TCARI/SAVI, and
TCARI/SAVI, and TCARI/OSAVI
TCARI/OSAVIduring during DOY
DOY 171–241
171–241 (20 (20 June—29
June—29 August
August 2017).
2017). Since
Since NDVI,
NDVI,
RDVI, SAVI, and OSAVI were closely related to the canopy structure; they increased
RDVI, SAVI, and OSAVI were closely related to the canopy structure; they increased with the growth with the growth
of maize
of maize and
and reached
reached thethe maximum
maximum on on DOY
DOY 207 207 (26
(26 July), decreased with
July), decreased with the
the gradual maturity of
gradual maturity of
maize. Previous studies showed that maize reaches full cover at a SAVI value of
maize. Previous studies showed that maize reaches full cover at a SAVI value of 0.64 [42,71]. SAVI of 0.64 [42,71]. SAVI of
each deficit irrigation region was greater than 0.61 on DOY 207, meaning maize
each deficit irrigation region was greater than 0.61 on DOY 207, meaning maize canopy had reached canopy had reached
effective coverage.
effective coverage. NDVI
NDVI had had the
the slowest
slowest decline
declineraterateandand the
the other
other three
three VIs
VIs were
were relatively
relatively close.
close.
Before the effective coverage of canopy (DOY 207), NDVI, RDVI, SAVI, and OSAVI
Before the effective coverage of canopy (DOY 207), NDVI, RDVI, SAVI, and OSAVI did not effectively did not effectively
distinguish deficit
distinguish deficit irrigation
irrigation treatments.
treatments. ForFor example,
example, During
During DOY
DOY 185 185 (the
(the start
start of
of deficit
deficit irrigation)
irrigation)
to 206,
to 206, the
the RDVI
RDVI values
values ofof TRT
TRT 1, 1, TRT
TRT 33 and
and TRT
TRT 55 were
were 0.76,
0.76, 0.74
0.74 and
and 0.73,
0.73, respectively,
respectively, with small
with small
difference no more than 0.03. The RDVI values of TRT 2 and TRT 4 were 0.69
difference no more than 0.03. The RDVI values of TRT 2 and TRT 4 were 0.69 and 0.67 with a small and 0.67 with a small
difference of
difference of 0.02.
0.02. There
There was
was clear
clear difference
difference between
between TRT TRT 33 and
and TRT
TRT 44 that
that had
had thethe same
same deficit
deficit
irrigation level (84%) with a difference value of 0.07. However, when the canopy
irrigation level (84%) with a difference value of 0.07. However, when the canopy reached effective reached effective
coverage, NDVI,
NDVI,RDVI,
RDVI,SAVI,SAVI, andandOSAVI
OSAVI started to show
started obvious
to show differences
obvious between
differences different
between deficit
different
irrigation
deficit treatments,
irrigation especially
treatments, in the maturation
especially in the maturationstage. stage.
For example, duringduring
For example, the maturation stage
the maturation
stage (DOY233~241), the RDVI values of TRT1 and TRT5 (100%), TRT2 and TRT3 (52%), and TRT4
(28%) were 0.54 and 0.54, 0.52 and 0.51, and 0.47, respectively. There was a clear numerical gradient
among different levels of deficit irrigation treatments.
In the case of TCARI, there were no significant changes during DOY 171–241. At the late
reproductive and maturation stages (DOY 231–241), different levels of deficit irrigation treatments
Remote Sens. 2019, 11, 605 12 of 24

Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 25

(DOY233~241), the RDVI values of TRT1 and TRT5 (100%), TRT2 and TRT3 (52%), and TRT4 (28%)
were coverage,
0.54 and 0.54,
they 0.52
couldandnot0.51, and 0.47,
effectively respectively.
distinguished There
water wasHowever,
stress. a clear numerical gradient
when canopy among
reached
different levels
effective of deficit
coverage, irrigation
they started totreatments.
distinguish deficit irrigation treatments gradually.

Figure 7. Changes
Figure 7. Changesof nine vegetation
of nine vegetationindices
indices(VIs):
(VIs): normalized difference
normalized difference vegetation
vegetation indexindex
(NDVI)(NDVI)
(a), renormalized difference
(a), renormalized vegetation
difference vegetationindex
index (RDVI) (b),soil-adjusted
(RDVI) (b), soil-adjusted vegetation
vegetation indexindex
(SAVI)(SAVI)
(c), (c),
optimization of soil-adjusted
optimization of soil-adjustedvegetation index(OSAVI)
vegetation index (OSAVI)(d),(d), transformed
transformed chlorophyll
chlorophyll absorption
absorption in
reflectance index
in reflectance index (TCARI)
(TCARI) (e),
(e), TCARI/NDVI
TCARI/NDVI (f), (f), TCARI/RDVI
TCARI/RDVI (g),(g),
TCARI/SAVI
TCARI/SAVI (h), (h),
and and
TCARI/OSAVI
TCARI/OSAVI (i), (i), during
during DOY
DOY 171–241(20
171–241 (20 June,
June, 20–29
20–29August)
August)forfor
allall
deficit irrigation
deficit treatments.
irrigation treatments.
DOY DOYis theisabbreviation
the abbreviation forfor
dayday
ofof year.The
year. Thered,
red, green,
green, blue,
blue,black
blackand
andpink lines
pink represent
lines the TRTs
represent the TRTs
1~5, respectively.
1~5, respectively.

3.4. Correlations between Vegetation Indices (VIs) and CWSI

In the case of TCARI, there were no significant changes during DOY 171–241. At the late
Tableand
reproductive 5 shows a description
maturation stagesof(DOY
correlations between
231–241), nine VIs
different and CWSI
levels (corresponding
of deficit scatter
irrigation treatments
plots are shown in Figures A2 and A3). When VIs and CWSI data from four days
could be distinguished gradually by TCARI. Similar to NDVI, RDVI, SAVI, and OSAVI, TCARI/NDVI,(DOY 223, 231, 240
and 241) inTCARI/SAVI
TCARI/RDVI, the reproductive
andand maturation stages
TCARI/OSAVI were
which used to perform
normalized linearRDVI,
by NDVI, regression
SAVIanalysis,
and OSAVI,
R2 were all less than 0.50. Among the four canopy structure VIs, NDVI had the best relationship (R2
changed with the growth of maize, reached minimum when the canopy reached effective coverage
= 0.34, p < 0.01). In the case of the five chlorophyll-related VIs, they all had a relatively better
(DOY 207), and then slowly increased. Before the canopy reached effective coverage, they could not
correlation (all R2 were greater than 0.35). Both TCARI/RDVI and TCARI/SAVI had significant
effectively distinguished
correlation with CWSIwater stress.
(p < 0.001) andHowever, when than
R2 were greater canopy
0.47.reached
When VIs effective
and CWSIcoverage, they
data from started
three
to distinguish deficit irrigation treatments gradually.
Remote Sens. 2019, 11, 605 13 of 24

3.4. Correlations between Vegetation Indices (VIs) and CWSI

Table 5 shows a description of correlations between nine VIs and CWSI (corresponding scatter
plots are shown in Figures A2 and A3). When VIs and CWSI data from four days (DOY 223, 231,
240 and 241) in the reproductive and maturation stages were used to perform linear regression analysis,
R2 were all less than 0.50. Among the four canopy structure VIs, NDVI had the best relationship
(R2 = 0.34, p < 0.01). In the case of the five chlorophyll-related VIs, they all had a relatively better
correlation (all R2 were greater than 0.35). Both TCARI/RDVI and TCARI/SAVI had significant
correlation with CWSI (p < 0.001) and R2 were greater than 0.47. When VIs and CWSI data from three
days (DOY 231, 240 and 241) in the late reproductive and maturation stages were used, the correlations
were significantly improved except for TCARI. All R2 were greater than 0.68 (p < 0.001). TCARI/RDVI
and TCARI/SAVI had the best relationship (both R2 ≥ 0.80, p < 0.001) Therefore, this paper choose
TCARI/RDVI and TCARI/SAVI to establish linear regression models between VI and CWSI. To rescale
the CWSI value between 0 and 1, the linear regression models were shown as Equations (14) and
(15), respectively.


 0 ( TCARI/RDV I ≤ 0.195)
CWSI − 1 = 2.41 ∗ ( TCARI/RDV I ) − 0.47 (0.195 < TCARI/RDV I < 0.609) (14)
(0.609 ≤ TCARI/RDV I )

 1


 0 ( TCARI/SAV I ≤ 0.182)
CWSI − 2 = 2.46 ∗ ( TCARI/SAV I ) − 0.45 (0.182 < TCARI/SAV I < 0.589) (15)
(0.589 ≤ TCARI/SAV I )

 1

Table 5. Coefficient of determination (R2 ), root mean square error (RMSE) from different predictions
with CWSI vs. nine VIs based on VIs and CWSI data from three (DOY 231, 240 and 241) or four (DOY
223, 231, 240 and 241) days.

Three Corresponding Data Four Corresponding Data

Vegetation Index
R2 (n = 15) RMSE R2 (n = 20) RMSE
Structural
NDVI 0.72 *** 0.046 0.34 ** 0.063
RDVI 0.81 *** 0.038 0.29 * 0.065
SAVI 0.81 *** 0.037 0.27 * 0.066
OSAVI 0.74 *** 0.044 0.29 * 0.065
Chlorophyll
TCARI 0.55 ** 0.058 0.35 ** 0.062
TCARI/NDVI 0.69 *** 0.048 0.42 ** 0.059
TCARI/RDVI 0.81 *** 0.037 0.47 *** 0.056
TCARI/SAVI 0.80 *** 0.039 0.50 *** 0.055
TCARI/OSAVI 0.77 *** 0.041 0.46 ** 0.057
* p < 0.05; ** p < 0.01; *** p < 0.001.

3.5. Maize Water Stress Maps Based on UAV Multispectral Remote-Sensing Imagery
Equations (14) and (15) were used to retrieve maize water stress maps based on UAV TCARI/RDVI
and TCARI/SAVI maps (Figure 8) for the reproductive (DOY 231) and maturation (DOY 240) stages.
Table 6 shows the comparison of CWSI values (at DOY 231 and DOY 240) using two different methods.
In the first method, CWSI values were calculated by on-site measurements within each sampling plot,
and in the second method, CWSI values, CWSI-1 and CWSI-2, retrieved from CWSI vs. VI regression
models (Equations (14) and (15)). The value of CWSI-1 and CWSI-2 was the average value of the entire
treatment region.
Remote Sens. 2019, 11, 605 14 of 24
Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 25

8. Maize
FigureFigure water
8. Maize stress
water maps
stress maps retrieved
retrievedbybyCWSI
CWSI vs. VI regression
vs. VI regressionmodels
models (Equations
(Equations (14) (14)
and and
(15)). (15)).
(a,b)(a)are
andwater
(b) arestress
watermaps retrieved
stress maps by Equation
retrieved by Equation (14)
(14)for
for the reproductive
the reproductive (DOY
(DOY 231)231)
and and
maturation
maturation (DOY (DOY
240)240) stages,
stages, respectively.
respectively. (c) and
(c,d) (d) are stress
are water water maps
stress maps retrieved
retrieved by Equation
by Equation (15) for
(15) for the reproductive (DOY 231) and maturation (DOY
the reproductive (DOY 231) and maturation (DOY 240) stages, respectively. 240) stages, respectively.

6. Comparison
TableTable of CWSI
6. Comparison of CWSIvalues (at(atDOY
values DOY231231 and
and DOY
DOY 240)240)using
usingtwotwo different
different methods.
methods. (1) (1)
CWSI values calculated by on-site measurements within each sampling plot, and
CWSI values calculated by on-site measurements within each sampling plot, and (2) CWSI values (2) CWSI values
retrieved
retrieved from from
CWSICWSI vs.regression
vs. VI VI regression models
models (Equations
(Equations (14)and
(14) and(15)).
(15)).CWSI-1
CWSI-1 and
and CWSI-2
CWSI-2 were
were the
the average value of the entire treatment
average value of the entire treatment region. region.
Reproductive (DOY 231) Maturation (DOY 240)
Treatment Reproductive (DOY
CWSI CWSI-1 CWSI-2231) CWSI Maturation CWSI-1 (DOY 240)CWSI-2
Treatment 0.06
TRT 1 0.09 0.09 0.02 0.11 0.11
CWSI CWSI-1 CWSI-2 CWSI CWSI-1 CWSI-2
TRT 2 0.04 0.08 0.08 0.20 0.18 0.18
TRT TRT
3 1 0.05 0.06 0.09 0.09 0.09 0.09 0.14 0.02 0.11
0.23 0.11
0.22
TRT TRT
4 2 0.19 0.04 0.19 0.08 0.18 0.08 0.26 0.20 0.40
0.18 0.39
0.18
TRT TRT
5 3 0.04 0.05 0.08 0.09 0.08 0.09 0.10 0.14 0.13
0.23 0.12
0.22
TRT 4 0.19 0.19 0.18 0.26 0.40 0.39
TRT 5 0.04 0.08 0.08 0.10 0.13 0.12
Remote Sens. 2019, 11, 605 15 of 24

4. Discussion
Previous research has reported that one of the drawbacks of CWSI empirical model was the
variability of NWSB [26,48,72,73]. Idso et al., (1981) reported that different crops had their own unique
NWSBs. Even for the same crop, NWSB will consequently change at different locations. For instance,
the slopes of maize NWSBs established at different locations vary from −1.10 to −3.77 ◦ C/kPa, and the
corresponding intercepts vary from 0.42 to 3.11 ◦ C. Most previous researchers ascribe this change to
the following factors:

(1) Climate variables [45,74–76], such as radiation and wind speed;

(2) Ta, RH measurement sitting [28,67,72], for example, Taghvaeian et al. [43] found that if Ta and
RH sensor was installed at full irrigation region, larger NWSB for sunflower would have been
assigned to deficit irrigation region;
(3) Tc monitoring instrument and its measure method [26,72];
(4) Specific hybrid differences [72,73], but previous studies have also suggested that the same NWSBs
can be used for similar varieties of maize [67] and wheat [77].

Many previous researchers have also observed that the coefficients of NWSB among annual crops
will change in different grown stages, such as maize, wheat, barley and sunflower. Idso et al. [69]
reported that the NWSBs of wheat and barley shifted significantly, as they move from the vegetative
to reproductive stage. The change of wheat and barley NWSBs was ascribed to the seed-heads that
absorbed the radiant energy, converted it into sensible heat, and provided a very effective aerodynamic
barrier against the transfer of sensible heat to the transpiring surfaces. Taghvaeian et al. [43] and
Nielsen [68] also observed significant differences in the coefficients of sunflower NWSBs during
the reproductive and maturation stages. The difference of sunflower NWSBs was attributed to the
decrease in sunflower transpiration and increase in head temperature after sunflower transitioned into
maturation stage. A similar observation was made by this study when to establish NWSBs for maize
during the reproductive and maturation stages.
When trying to explain the difference of NWSBs among annual crops in different grown stages,
previous researchers focused on the change of crop physiological characteristics, but ignored the change
of microclimate, such as Ta and RH. In this study, compared to the reproductive stage, the average
daily maximum temperature in the maturation stage decreased from 31.31 ◦ C to 25.46 ◦ C and average
daily minimum relative humidity increased from 29.78% to 33.23%, with a decrease of 5.85 ◦ C and
an increase of 3.45%, respectively (Table 2). The average daily maximum temperature and average
daily minimum relative humidity had a similar changing trend in the study of Taghvaeian et al. [43].
Figure 9a shows the changing trend of VPD along with Ta and RH. The increase of RH and decrease of
Ta will lead to a decrease of VPD. This result will lead to a decrease in transpiration rate of crops under
well-watered conditions, which means the difference between Tc and Ta will lower. In other words,
the decrease in VPD will lead to an increase in NWSB. In this study, compared to the reproductive
stage, VPD decreased in the maturation stage with average value of 0.8 kPa, and NWSB increased with
average value of 2.1 ◦ C.
In general terms, plants under well irrigation condition have specific non-water stress lines
that corresponding to average CWSI values. Previous studies about the use of the CWSI as a water
stress indicator in olive and vineyards orchards have indicated that average CWSI in a well irrigation
condition was less than 0.5 [4,16,40,78]. For maize, in this study, average CWSI values were no
more than 0.12 in the reproductive stage, and a much smaller average CWSI, a value of 0.03, in the
maturation stage. An interesting result could be found that average CWSI value for maize under full
irrigation condition, the non-water stress line, was much smaller, compared to the reproductive stage.
The same average CWSI calculated by on- site measurements represented more serious water stress
in the maturation stage than the reproductive stage. For example, the average CWSI of TRT 2 in the
reproductive stage was 0.12, while the average CWSI of TRT 3 in the maturation stage was also 0.12,
with 45% water applied difference. The average CWSI of TRT 4 in the reproductive stage was 0.27,
 [69] reported that the NWSBs of wheat and barley shifted significantly, as they move from the
vegetative to reproductive stage. The change of wheat and barley NWSBs was ascribed to the seed-
heads that absorbed the radiant energy, converted it into sensible heat, and provided a very effective
aerodynamic barrier against the transfer of sensible heat to the transpiring surfaces. Taghvaeian et al.
[43] and Nielsen [68] also observed significant differences in the coefficients of sunflower NWSBs
Remote Sens. 2019, 11, 605 16 of 24
during the reproductive and maturation stages. The difference of sunflower NWSBs was attributed
to the decrease in sunflower transpiration and increase in head temperature after sunflower
transitioned into maturation stage. A similar observation was made by this study when to establish
while the average CWSI of TRT 4 in the maturation stage was 0.32. 43% water applied difference
NWSBs for maize during the reproductive and maturation stages.
corresponded When totrying
a 0.05todifference
explain theof CWSI value.
difference of NWSBsSimilar
amongphenomenon
annual crops had been observed
in different by Han
grown stages,
Ming et al. [26]researchers
previous for maize in Greeley,
focused CO,change
on the USA. In of the
cropreproductive
physiological stage, average but
characteristics, CWSI values
ignored thewere
no more thanof0.16,
change and 0.5 insuch
microclimate, the maturation
as Ta and RH.stage. This
In this phenomenon
study, compared to may
thebe explained stage,
reproductive in twotheways:
average daily maximum temperature in the maturation stage decreased from 31.31 °C to 25.46 °C and
(1) The change of maize physiological characteristics and microclimate from the reproductive to the
average daily minimum relative humidity increased from 29.78% to 33.23%, with a decrease of 5.85
maturation stage causes significant difference in NWSBs and insignificant difference in NTBs.
°C and an increase of 3.45%, respectively (Table 2). The average daily maximum temperature and
(2) The cropdaily
average and minimum
soil characteristics in the field
relative humidity had ahave significant
similar changingvariability
trend in theinstudy
eachofdeficit irrigation
Taghvaeian
treatment
et al. [43]. Figure 9a shows the changing trend of VPD along with Ta and RH. The increase of RH and 8).
regions. One sample plot could not represent the entire treatment region well (Figure
Santesteban
decrease of Ta et will
al. [79]
leadhad
to aalso foundofthat
decrease VPD.anThis
18-year-old
result willvineyard
lead to aunder
decreasesame irrigation condition
in transpiration rate
of crops
and under
without well-watered
significant petsconditions, which
and diseases hadmeans the difference
significant between
variability Tc and
in water Ta will lower.
status.
In other words, the decrease in VPD will lead to an increase in NWSB. In this study, compared to the

(a) (b)

Figure 9. Vapor pressure deficit (VPD) changing trend along with Ta and RH (a); Changes of NWSB
between the reproductive and maturation stages (b).

Previous studies in the vineyard have indicated that multispectral VIs acquired from
remote-sensing images have a significantly high correlation with water stress indicators, revealing
the potential to assess water stress and its variability by using multispectral VIs. Espinoza et al. [80]
found a high correlation (r = 0.65, p < 0.01) between GNDVI and stomatal conductance. Baluja et
al. [81] found a R2 value of 0.68 (n = 10) from the linear regression between NDVI and stem water
potential. Zarco-Tejada et al. [33] found a R2 value of 0.77 (n = 9) between PRInorm and CWSI. However,
the relationship between UAV multispectral VIs and CWSI is not well documented for maize.
In this study, a higher significantly correlation between VIs acquired from UAV multispectral
images and CWSI was found. In particular, compared to CWSI, the R2 values of TCARI/RDVI
and TCARI/SAVI were 0.81 and 0.80 (n = 15) in the late reproductive and maturation stages.
A similar observation was made by Baluja et al. [81] in a commercial vineyard, with the R2 values of
TCARI/OSAVI were 0.58 and 0.84 (n = 10) compared to stem water potential and stomatal conductance,
respectively. The results showed that TCARI/OSAVI were probably reflecting the result of cumulative
water deficits, and hence a long-term response. In this study similar results were found, namely,
when using UAV multispectral VIs to estimate maize water stress, their effects was influenced by the
extended period of water stress. With the prolongation of water stress, UAV multispectral VIs could
better distinguish different levels of deficit irrigation treatments (Figure 7). A similar phenomenon
was also found by Espinoza et al. [80], Li et al. [82] and Zulini et al. [83]. Perhaps, this could explain
that the R2 values of VIs-CWSI were higher in the late reproductive and maturation stages (DOY 231,
240 and 241) than in the reproductive and maturation stages (DOY 223, 231, 240 and 241) (Table 5).
However, Zarco-Tejada et al. [84] found that TCARI/OSAVI was less correlated with leaf stomatal
conductance and water potential in Citrus orchards, with the R2 values of 0.45 and 0.51. The reason for
the different observations might be that TCARI/OSAVI are more related to plant vigor than the plant
dynamic physiological status, with the result that it might correlate well in crops where the biomass
Remote Sens. 2019, 11, 605 17 of 24

proportionally increases in parallel to photosynthesis, however, in the case of evergreen crop species,
results were not so satisfactory [85].
TCARI/OSAVI, similar to TCARI/RDVI and TCARI/SAVI in this study, was established by
Haboudane et al. [60] to make accurate predictions of maize chlorophyll content from hyperspectral
remote-sensing imagery in Canada. The results showed that TCARI/OSAVI was both very sensitive to
chlorophyll content variations and very resistant to the variations of leaf area index (LAI). The results of
Zarco-Tejada et al. [86] in vineyard showed that TCARI/OSVAI was also correlated with xanthophyll
and the carotenoid content with R2 value of 0.51 (p < 0.01). Under water stress condition, transpiration
rate of crop may reduce, resulting a reduction of cooling effect [87,88]. To avoid damage to the
tissues, crops need to reduce the absorption of light and dissipate the extra energy through the
reduction of chlorophyll b and the interconversion of the xanthophyll cycle pigments, respectively [89].
Therefore, the correlations of TCARI/OSAVI with chlorophyll and xanthophyll content may be the
crop physiological basis of why TCARI/OSAVI family of indices could monitor crop water stress
status. Future, related research is still needed.
From the CWSI maps retrieved by TCARI/RDVI and TCARI/SAVI, we could find that the water
stress status of maize in each irrigation treatment had spatial variations. Compared to CWSI calculated
by on-site measurements, CWSI based on VI-CWSI regression models could better reflect water stress
conditions of maize in the field scale. Taking CWSI-2 as an example, the average CWSI-2 could
reflect TRT 4 (69%) and TRT 1(100%) in the reproductive stage, with the values of 0.18 and 0.09,
respectively; and could reflect TRT4 (28%), TRT2 and TRT3 (52%), and TRT1 (100%) in the maturation
stage, with the value of 0.39, 0.18 and 0.22, and 0.11, respectively. Compared with CWSI calculated by
on-site measurements, CWSI-2 was more reasonable (Table 6).

5. Conclusions
As one of the most important remote-sensing platforms, an unmanned aerial vehicle (UAV)
system could obtain key information (e.g., water stress status) about a crop at high temporal frequency
(1–3 days revisit) and high spatial resolution (centimeter scale). This study investigated whether
a UAV-based multispectral remote sensing system could map the water stress status of maize at a
farm scale, as a supplement to ground-based empirical crop water stress index (CWSI) model. Since,
the local weather condition and management practices can greatly influence the applicability of the
empirical CWSI model, changes of the non-water stress baseline (NWSB) and the non-transpiring
baseline (NTB) were also analyzed. Our results confirmed that NWSB had a significant difference in
the reproductive and maturation stages, however NTB did not change significantly. The shift of NWSB
from the reproductive to maturation stage is not only related to the change of maize characteristics
(e.g., canopy structure and transpiration rate), but also related to the change of microclimate, such
as air temperature and relative humidity. The ratio of transformed chlorophyll absorption in the
reflectance index (TCARI) and renormalized difference vegetation index (RDVI), and the TCARI
and soil-adjusted vegetation index (SAVI) had high correlations with CWSI. R2 values were 0.47 and
0.50 for TCARI/RDVI at the reproductive and maturation stages, respectively; and 0.81 and 0.80 for
TCARI/SAVI at the late reproductive and maturation stages, respectively. This study demonstrates the
feasibility of mapping maize crop water status and monitoring its spatial variability at a farm scale by
using UAV-based multispectral VI-CWSI regression models established in this study. As a future work,
research into the influence of other meteorological parameters (e.g., extreme illumination) on the water
stress monitoring effect of CWSI will be conducted. What changes of crop physiological characteristics
(e.g., the changes of chlorophyll and xanthophyll caused by photoprotection) enable the vegetation
index (e.g., TCARI/SAVI) to be used to monitor the water stress status of maize will also be studied.

Author Contributions: L.Z. and W.H. conceived and designed the experiments; L.Z. and Y.N. analyzed the data;
L.Z. and H.Z. discussed and drafted the manuscript. H.Z. and W.H. revised the manuscript and edited English
language. All authors read and approved the final version.
Remote Sens. 2019, 11, 605 18 of 24

Funding: This study was supported by the National Key R & D plan from the MOST of China (2017YFC0403203),
the Synergetic Innovation of Industry-University-Research Cooperation Project plan from Yangling (2018CXY-23),
the 111 Project (No.B12007) and the Key Discipline Construction Project of Northwest Agriculture and Forestry
University (2017-C03).
Acknowledgments: We are grateful to Guomin Shao, Yi Wang and Guang Li for data collection.
Remote Sens. 2019, 11, x FOR PEER REVIEW 18 of 25
Conflicts of Interest: No conflict of interest.
Conflicts of Interest: No conflict of interest.
Appendix A
Appendix A

Figure A1. TheA1.

Figure changing curvescurves
The changing of meteorological datadata
of meteorological during
duringthethe
reproductive
reproductiveand andmaturation
maturation stages
in 2017. stages
(a): airintemperature (Ta); (b) relative
2017. (a): air temperature humidity
(Ta); (b) (RH); (c) (RH);
relative humidity net solar radiation
(c) net (Rn); (d)
solar radiation weed
(Rn); (d) speed.
weedabbreviation
DOY is the speed. DOY isfor the day
abbreviation
of year.for dayblue,
The of year.
red The blue,
and red and
black black
solid solid
lines lines represent
represent the
the maximum,
minimum maximum,
and mean minimum
values,and mean values,The
respectively. respectively. The black
black dotted line dotted line is boundary
is boundary betweenbetween the
the reproductive
reproductive and maturation stages. The blue dotted lines indicate dates when data was taken for the
and maturation stages. The blue dotted lines indicate dates when data was taken for the calculation of
calculation of non-water-stressed and non-transpiring baselines of CWSI.
non-water-stressed and non-transpiring baselines of CWSI.
Remote Sens. 2019, 11, 605 19 of 24
Remote Sens. 2019, 11, x FOR PEER REVIEW 19 of 25

Figure Correlations
A2.A2.
Figure Correlationsbetween
betweenunmanned
unmannedaerial
aerial vehicle (UAV)-basedmultispectral
vehicle (UAV)-based multispectralVIsVIsand
and CWSI
CWSI
based onon
three (DOY 231, 2 represents coefficient of determination.
based three (DOY 231,240
240and
and241)
241)corresponding
corresponding data. RR2 represents coefficient of determination.
p represents Pearson
p represents correlation
Pearson correlationcoefficient.
coefficient.
Remote Sens. 2019, 11, 605 20 of 24

Remote Sens. 2019, 11, x FOR PEER REVIEW 20 of 25

Figure Correlations
A3.A3.
Figure Correlations between
betweenUAV-based
UAV-based multispectral
multispectralVIs
VIsand CWSI
and CWSIbased
basedononour (DOY
our (DOY223, 231,
223,
2
240231,
and240
241)
andcorresponding data.data.
241) corresponding R represents coefficient
R2 represents determination. pprepresents
coefficientofofdetermination. representsPearson
Pearson
correlation coefficient.
correlation coefficient.

References
References
1. 1 Lopez, J.R.;
Lopez, Winter,
J.R.; J.M.;
Winter, Elliott,
J.M.; Elliott,J.;J.;Ruane,
Ruane,A.C.;
A.C.;Porter,
Porter, C.; Hoogenboom,
Hoogenboom,G. G.Integrating
Integratinggrowth
growth stage
stage deficit
deficit
irrigation into into
irrigation a process
a based
process cropbased Agric. model.
model.crop For. Meteorol.
Agric.2017, 243,Meteorol.
For. 84–92. [CrossRef]
2017, 243, 84–92,
2. Romero-Trigueros, C.; Nortes, P.A.; Alarcon, J.J.; Hunink, J.E.; Parra, M.; Contreras, S.; Droogers, P.; Nicolas, E.
doi:10.1016/j.agrformet.2017.05.001.
2 Effects
Romero-Trigueros,
of saline reclaimedC.; Nortes,
waters and P.A.; Alarcon,
deficit J.J.; Hunink,
irrigation on CitrusJ.E.; Parra, M.;
physiology Contreras,
assessed S.; Droogers,
by UAV P.;
remote sensing.
Nicolas,
Agric. WaterE.Manag.
Effects 2017,
of saline
183,reclaimed waters and deficit irrigation on Citrus physiology assessed by UAV
60–69. [CrossRef]
3. remote
Zhao, T.B.;sensing. Agric.
Stark, B.; Chen,Water Manag.
Y.Q.; 2017, Doll,
Ray, A.L.; 183, 60–69, doi:10.1016/j.agwat.2016.09.014.
D. A Detailed Field Study of Direct Correlations Between
3 Ground
Zhao,Truth
T.B.; Stark, B.; Chen, Y.Q.; Ray, A.L.; Doll,
Crop Water Stress and Normalized Difference D. A Detailed Field Study
Vegetation Indexof(NDVI)
Direct Correlations Between
from Small Unmanned
Ground Truth Crop Water Stress and Normalized Difference Vegetation Index
Aerial System (sUAS). In Proceedings of the 2015 International Conference on Unmanned Aircraft Systems (NDVI) from Small
Unmanned
(ICUAS), Denver,Aerial
CO,System (sUAS).
USA, 9–12 JuneIn2015;
Proceedings
pp. 520–525.of the 2015 International Conference on Unmanned
4. Sepulveda-Reyes, D.; Ingram, B.; Bardeen, M.; Zuniga, M.;2015;
Aircraft Systems (ICUAS), Denver, CO, USA, 9–12 June pp. 520–525.
Ortega-Farias, S.; Poblete-Echeverria, C. Selecting
4 Sepulveda-Reyes, D.; Ingram, B.; Bardeen, M.; Zuniga, M.; Ortega-Farias, S.; Poblete-Echeverria, C.
Canopy Zones and Thresholding Approaches to Assess Grapevine Water Status by Using Aerial and
Selecting Canopy Zones and Thresholding Approaches to Assess Grapevine Water Status by Using Aerial
Ground-Based Thermal Imaging. Remote Sens. 2016, 8, 822. [CrossRef]
and Ground-Based Thermal Imaging. Remote Sens. 2016, 8, 822, doi:10.3390/rs8100822.
5. Cui, N.; Du, T.; Kang, S.; Li, F.; Zhang, J.; Wang, M.; Li, Z. Regulated deficit irrigation improved fruit quality
and water use efficiency of pear-jujube trees. Agric. Water Manag. 2008, 95, 489–497. [CrossRef]
Remote Sens. 2019, 11, 605 21 of 24

6. Geerts, S.; Raes, D. Deficit irrigation as an on-farm strategy to maximize crop water productivity in dry areas.
Agric. Water Manag. 2009, 96, 1275–1284. [CrossRef]
7. Han, M.; Zhang, H.; DeJonge, K.C.; Comas, L.H.; Trout, T.J. Estimating maize water stress by standard
deviation of canopy temperature in thermal imagery. Agric. Water Manag. 2016, 177, 400–409. [CrossRef]
8. Ihuoma, S.O.; Madramootoo, C.A. Recent advances in crop water stress detection. Comput. Electron. Agric.
2017, 141, 267–275. [CrossRef]
9. Li, L.; Nielsen, D.C.; Yu, Q.; Ma, L.; Ahuja, L.R. Evaluating the Crop Water Stress Index and its correlation
with latent heat and CO2 fluxes over winter wheat and maize in the North China plain. Agric. Water Manag.
2010, 97, 1146–1155. [CrossRef]
10. Campbell, G.S.; Campbell, M.D.J.A.i.I. Irrigation Scheduling Using Soil Moisture Measurements: Theory
and Practice. Adv. Irrig. 1982, 1, 25–42. [CrossRef]
11. Hazaymeh, K.; Hassan, Q.K. Remote sensing of agricultural drought monitoring: A state of art review.
Aims Environ. Sci. 2016, 3, 604–630. [CrossRef]
12. Calera, A.; Campos, I.; Osann, A.; D’Urso, G.; Menenti, M. Remote Sensing for Crop Water Management:
From ET Modelling to Services for the End Users. Sensors 2017, 17, 1104. [CrossRef]
13. Du, L.; Tian, Q.; Yu, T.; Meng, Q.; Jancso, T.; Udvardy, P.; Huang, Y. A comprehensive drought monitoring
method integrating MODIS and TRMM data. Int. J. Appl. Earth Obs. Geoinf. 2013, 23, 245–253. [CrossRef]
14. Helman, D.; Bahat, I.; Netzer, Y.; Ben-Gal, A.; Alchanatis, V.; Peeters, A.; Cohen, Y. Using Time Series
of High-Resolution Planet Satellite Images to Monitor Grapevine Stem Water Potential in Commercial
Vineyards. Remote Sens. 2018, 10, 1615. [CrossRef]
15. Zhang, H.; Han, M.; Chavez, J.L.; Lan, Y. Improvement in estimation of soil water deficit by integrating
airborne imagery data into a soil water balance model. Int. J. Agric. Biol. Eng. 2017, 10, 37–46. [CrossRef]
16. Bellvert, J.; Marsal, J.; Girona, J.; Zarco-Tejada, P.J. Seasonal evolution of crop water stress index in grapevine
varieties determined with high-resolution remote sensing thermal imagery. Irrig. Sci. 2015, 33, 81–93.
[CrossRef]
17. Gerhards, M.; Schlerf, M.; Rascher, U.; Udelhoven, T.; Juszczak, R.; Alberti, G.; Miglietta, F.; Inoue, Y. Analysis
of Airborne Optical and Thermal Imagery for Detection of Water Stress Symptoms. Remote Sens. 2018, 10,
1139. [CrossRef]
18. Bai, H.; Purcell, L.C. Aerial canopy temperature differences between fast- and slow-wilting soya bean
genotypes. J. Agron. Crop Sci. 2018, 204, 243–251. [CrossRef]
19. Herwitz, S.R.; Johnson, L.F.; Dunagan, S.E.; Higgins, R.G.; Sullivan, D.V.; Zheng, J.; Lobitz, B.M.; Leung, J.G.;
Gallmeyer, B.A.; Aoyagi, M.; et al. Imaging from an unmanned aerial vehicle: Agricultural surveillance and
decision support. Comput. Electron. Agric. 2004, 44, 49–61. [CrossRef]
20. Han, W.; Zhang, L.; Zhang, H.; Shi, Z.; Yuan, M.; Wang, Z. Extraction Method of Sublateral Canal Distribution
Information Based on UAV Remote Sensing. Nongye Jixie Xuebao/Trans. Chin. Soc. Agric. Mach. 2017, 48,
210–219. [CrossRef]
21. Han, W.; Guang, L.I.; Yuan, M.; Zhang, L.; Shi, Z. Extraction Method of Maize Planting Information Based
on UAV Remote Sensing Techonology. Nongye Jixie Xuebao/Trans. Chin. Soc. Agric. Mach. 2017, 48, 139–147.
[CrossRef]
22. Zhao, C. Advances of Research and Application in Remote Sensing for Agriculture. Nongye Jixie Xuebao/Trans.
Chin. Soc. Agric. Mach. 2014, 45, 277–293. [CrossRef]
23. Inoue, Y.; Kimball, B.A.; Jackson, R.D.; Jr, P.J.P.; Reginato, R.J. Remote estimation of leaf transpiration rate
and stomatal resistance based on infrared thermometry. Agric. For. Meteorol. 1990, 51, 21–33. [CrossRef]
24. Zhang, Z.; Bian, J.; Han, W.; Fu, Q.; Chen, S.; Cui, T. Cotton moisture stress diagnosis based on canopy
temperature characteristics calculated from UAV thermal infrared image. Nongye Gongcheng Xuebao/Trans.
Chin. Soc. Agric. Eng. 2018, 34, 77–84. [CrossRef]
25. Martínez, J.; Egea, G.; Agüera, J.; Pérez-Ruiz, M. A cost-effective canopy temperature measurement system
for precision agriculture: A case study on sugar beet. Precis. Agric. 2016, 18, 95–110. [CrossRef]
26. Han, M.; Zhang, H.H.; DeJonge, K.C.; Comas, L.H.; Gleason, S. Comparison of three crop water stress index
models with sap flow measurements in maize. Agric. Water Manag. 2018, 203, 366–375. [CrossRef]
27. Zia, S.; Romano, G.; Spreer, W.; Sanchez, C.; Cairns, J.; Araus, J.L.; Muller, J. Infrared Thermal Imaging as a
Rapid Tool for Identifying Water-Stress Tolerant Maize Genotypes of Different Phenology. J. Agron. Crop Sci.
2013, 199, 75–84. [CrossRef]
Remote Sens. 2019, 11, 605 22 of 24

28. Zhang, L.; Niu, Y.; Han, W.; Liu, Z. Establishing Method of Crop Water Stress Index Empirical Model of Field
Maize. Nongye Jixie Xuebao/Trans. Chin. Soc. Agric. Mach. 2018, 49, 233–239. [CrossRef]
29. Irmak, S.; Haman, D.Z.; Bastug, R. Determination of crop water stress index for irrigation timing and yield
estimation of corn. Agron. J. 2000, 92, 1221–1227. [CrossRef]
30. Zhang, Z.; Bian, J.; Han, W.; Fu, Q.; Chen, S.; Cui, T. Diagnosis of Cotton Water Stress Using Unmanned
Aerial Vehicle Thermal Infrared Remote Sensing after Removing Soil. Nongye Jixie Xuebao/Trans. Chin. Soc.
Agric. Mach. 2018, 49, 250–260. [CrossRef]
31. Cohen, Y.; Alchanatis, V.; Sela, E.; Saranga, Y.; Cohen, S.; Meron, M.; Bosak, A.; Tsipris, J.; Ostrovsky, V.;
Orolov, V.; et al. Crop water status estimation using thermography: Multi-year model development using
ground-based thermal images. Precis. Agric. 2015, 16, 311–329. [CrossRef]
32. Pou, A.; Diago, M.P.; Medrano, H.; Baluja, J.; Tardaguila, J. Validation of thermal indices for water status
identification in grapevine. Agric. Water Manag. 2014, 134, 60–72. [CrossRef]
33. Zarco-Tejada, P.J.; Gonzalez-Dugo, V.; Williams, L.E.; Suarez, L.; Berni, J.A.J.; Goldhamer, D.; Fereres, E.
A PRI-based water stress index combining structural and chlorophyll effects: Assessment using diurnal
narrow-band airborne imagery and the CWSI thermal index. Remote Sens. Environ. 2013, 138, 38–50.
[CrossRef]
34. Agam, N.; Cohen, Y.; Alchanatis, V.; Ben-Gal, A. How sensitive is the CWSI to changes in solar radiation?
Int. J. Remote Sens. 2013, 34, 6109–6120. [CrossRef]
35. Bellvert, J.; Zarco-Tejada, P.J.; Girona, J.; González-Dugo, V.; Fereres, E. A tool for detecting crop water status
using airborne high-resolution thermal imagery. WIT Trans. Ecol. Environ. 2014, 185, 25–32.
36. Wang, D.; Gartung, J. Infrared canopy temperature of early-ripening peach trees under postharvest deficit
irrigation. Agric. Water Manag. 2010, 97, 1787–1794. [CrossRef]
37. Paltineanu, C.; Septar, L.; Moale, C. Crop Water Stress in Peach Orchards and Relationships with Soil
Moisture Content in a Chernozem of Dobrogea. J. Irrig. Drain. Eng. 2013, 139, 20–25. [CrossRef]
38. Agam, N.; Cohen, Y.; Berni, J.A.J.; Alchanatis, V.; Kool, D.; Dag, A.; Yermiyahu, U.; Ben-Gal, A. An insight to
the performance of crop water stress index for olive trees. Agric. Water Manag. 2013, 118, 79–86. [CrossRef]
39. Gençǧlan, C.; Yazar, A. Determination of Crop Water Stress İndex (CWSI) and irrigation timing by utilizing
infrared thermometer values on the first corn grown under çukurova conditions. Am. J. Surg. 2014, 3,
342–345. [CrossRef]
40. Berni, J.A.J.; Zarco-Tejada, P.J.; Sepulcre-Canto, G.; Fereres, E.; Villalobos, F. Mapping canopy conductance
and CWSI in olive orchards using high resolution thermal remote sensing imagery. Remote Sens. Environ.
2009, 113, 2380–2388. [CrossRef]
41. DeJonge, K.C.; Taghvaeian, S.; Trout, T.J.; Comas, L.H. Comparison of canopy temperature-based water
stress indices for maize. Agric. Water Manag. 2015, 156, 51–62. [CrossRef]
42. Taghvaeian, S.; Chavez, J.L.; Hansen, N.C. Infrared Thermometry to Estimate Crop Water Stress Index and
Water Use of Irrigated Maize in Northeastern Colorado. Remote Sens. 2012, 4, 3619–3637. [CrossRef]
43. Taghvaeian, S.; Comas, L.; DeJonge, K.C.; Trout, T.J. Conventional and simplified canopy temperature indices
predict water stress in sunflower. Agric. Water Manag. 2014, 144, 69–80. [CrossRef]
44. Wang, L.M.; Qiu, G.Y.; Zhang, X.Y.; Chen, S.Y. Application of a new method to evaluate crop water stress
index. Irrig. Sci. 2005, 24, 49–54. [CrossRef]
45. Gonzalez-Dugo, V.; Zarco-Tejada, P.J.; Fereres, E. Applicability and limitations of using the crop water stress
index as an indicator of water deficits in citrus orchards. Agric. For. Meteorol. 2014, 198, 94–104. [CrossRef]
46. Ballester, C.; Jimenez-Bello, M.A.; Castel, J.R.; Intrigliolo, D.S. Usefulness of thermography for plant water
stress detection in citrus and persimmon trees. Agric. For. Meteorol. 2013, 168, 120–129. [CrossRef]
47. Gonzalez-Dugo, V.; Goldhamer, D.; Zarco-Tejada, P.J.; Fereres, E. Improving the precision of irrigation in a
pistachio farm using an unmanned airborne thermal system. Irrig. Sci. 2015, 33, 43–52. [CrossRef]
48. Idso, S.B.; Jackson, R.D.; Pinter, P.J., Jr.; Reginato, R.J.; Hatfield, J.L. Normalizing the stress-degree-day
parameter for environmental variability. Agric. Meteorol. 1981, 24, 45–55. [CrossRef]
49. Jackson, R.D.; Idso, S.B.; Reginato, R.J.; Pinter, P.J. Canopy temperature as a crop water stress indicator. Water
Resour. Res. 1981, 17, 1133–1138. [CrossRef]
50. Rodriguez-Perez, J.R.; Riano, D.; Carlisle, E.; Ustin, S.; Smart, D.R. Evaluation of hyperspectral reflectance
indexes to detect grapevine water status in vineyards. Am. J. Enol. Vitic. 2007, 58, 302–317.
Remote Sens. 2019, 11, 605 23 of 24

51. Park, S.; Ryu, D.; Fuentes, S.; Chung, H.; Hernández-Montes, E.; O’Connell, M. Adaptive Estimation of Crop
Water Stress in Nectarine and Peach Orchards Using High-Resolution Imagery from an Unmanned Aerial
Vehicle (UAV). Remote Sens. 2017, 9, 828. [CrossRef]
52. Moran, M.S.; Clarke, T.R.; Inoue, Y.; Vidal, A. Estimating Crop Water-Deficit Using the Relation Between
Surface-Air Temperature and Spectral Vegetation Index. Remote Sens. Environ. 1994, 49, 246–263. [CrossRef]
53. Wang, P.; Luo, X.; Zhou, Z.; Zang, Y.; Hu, L. Key technology for remote sensing information acquisitionbased
on micro UAV. Nongye Gongcheng Xuebao/Trans. Chin. Soc. Agric. Eng. 2014, 30, 1–12. [CrossRef]
54. Xiang, H.; Tian, L. Development of a low-cost agricultural remote sensing system based on an autonomous
unmanned aerial vehicle (UAV). Biosyst. Eng. 2011, 108, 174–190. [CrossRef]
55. Bellvert, J.; Zarco-Tejada, P.J.; Girona, J.; Fereres, E. Mapping crop water stress index in a ’Pinot-noir’ vineyard:
Comparing ground measurements with thermal remote sensing imagery from an unmanned aerial vehicle.
Precis. Agric. 2014, 15, 361–376. [CrossRef]
56. Ribeiro-Gomes, K.; Hernandez-Lopez, D.; Ortega, J.F.; Ballesteros, R.; Poblete, T.; Moreno, M.A. Uncooled
Thermal Camera Calibration and Optimization of the Photogrammetry Process for UAV Applications in
Agriculture. Sensors 2017, 17, 2173. [CrossRef]
57. Ni, G. Vegetation Index and Its Advances. Arid Meteorol. 2003, 21, 71–75.
58. Rouse, J.W., Jr.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring Vegetation Systems in the Great Plains
with Erts. Nasa Spec. Publ. 1973, 351, 309.
59. Huete, A.R. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ. 1988, 25, 295–309. [CrossRef]
60. Driss, H.; John, R.M.; Nicolas, T.; Pablo, J.Z.-T.; Louise, D. Integrated narrow-band vegetation indices for
prediction of crop chlorophyll content for application to precision agriculture. Remote Sens. Environ. 2002, 81,
416–426. [CrossRef]
61. Lv, Y.; Li, B. Soil Science; China Agriculture Press: Beijing, China, 2006.
62. Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop Evapotranspiration—Guidelines for Computing Crop Water
Requirements—FAO Irrigation and Drainage Paper 56; FAO—Food and Agriculture Organization of the United
Nations: Rome, Italy, 1998.
63. Li, Y.; Yan, H.; Cai, D.; Gu, T.; Sui, R.; Chen, D. Evaluating the water application uniformity of center pivot
irrigation systems in Northern China. Int. Agric. Eng. J. 2018. under review.
64. Heermann, D.F.H.; Hein, P.R. Performance characteristics of self-propelled center-pivot sprinkler irrigation
system. Trans. ASAE 1968, 11, 11–15.
65. Yang, G.; Li, C.; Yu, H.; Xu, B.; Feng, H.; Gao, L.; Zhu, D. UAV based multi-load remote sensing technologies
for wheat breeding information acquirement. Nongye Gongcheng Xuebao/Trans. Chin. Soc. Agric. Eng. 2015,
31, 184–190. [CrossRef]
66. Walter, I.A.; Allen, R.G.; Elliott, R.; Jensen, M.E.; Itenfisu, D.; Mecham, B.; Howell, T.A.; Snyder, R.; Brown, P.;
Echings, S. ASCE’s standardized reference evapotranspiration equation. In Proceedings of the Watershed
Management and Operations Management 2000, Fort Collins, CO, USA, 20–24 June 2000; pp. 1–11.
67. Taghvaeian, S.; Chávez, J.L.; Bausch, W.C.; DeJonge, K.C.; Trout, T.J. Minimizing instrumentation requirement
for estimating crop water stress index and transpiration of maize. Irrig. Sci. 2014, 32, 53–65. [CrossRef]
68. Nielsen, D.C. Non water-stressed baselines for sunflowers. Agric. Water Manag. 1994, 26, 265–276. [CrossRef]
69. Idso, S.B. Non-water-stressed baselines: A key to measuring and interpreting plant water stress. Agric.
Meteorol. 1982, 27, 59–70. [CrossRef]
70. Xiao, C.; Xu, L.; Yuan, G.; Wang, W.; Yi, L. Crop water stress index model for monitoring summer maize
water stress based on canopy surface temperature. Trans. Chin. Soc. Agric. Eng. 2005, 21, 22–24. [CrossRef]
71. Bausch, W.C. Soil background effects on reflectance-based crop coefficients for corn. Remote Sens. Environ.
1993, 46, 213–222. [CrossRef]
72. Yazar, A.; Howell, T.A.; Dusek, D.A.; Copeland, K.S. Evaluation of crop water stress index for LEPA irrigated
corn. Irrig. Sci. 1999, 18, 171–180. [CrossRef]
73. Gardner, B.R.; Nielsen, D.C.; Shock, C.C. Infrared thermometry and the crop water stress index. II. Sampling
procedures and interpretation. J. Prod. Agric. 1993, 5, 466–475. [CrossRef]
74. Payero, J.O.; Irmak, S. Variable upper and lower crop water stress index baselines for corn and soybean.
Irrig. Sci. 2006, 25, 21–32. [CrossRef]
75. Payero, J.O.; Neale, C.M.U.; Wright, J.L. Non-water-stressed baselines for calculating Crop Water Stress
Index (CWSI) for alfalfa and tall fescue grass. Trans. ASAE 2005, 48, 653–661. [CrossRef]
Remote Sens. 2019, 11, 605 24 of 24

76. Zolnier, S.; Gates, R.S.; Anderson, R.G.; Nokes, S.E.; Duncan, G.A. Non-water-stressed baseline as a tool
for dynamic control of a misting system for propagation of poinsettias. Trans. ASAE 2001, 44, 137–147.
[CrossRef]
77. Idso, S.B.; Reginato, R.J.; Clawson, K.L.; Anderson, M.G. On the stability of non-water-stressed baselines.
Agric. For. Meteorol. 1984, 32, 177–182. [CrossRef]
78. Ben-Gal, A.; Agam, N.; Alchanatis, V.; Cohen, Y.; Yermiyahu, U.; Zipori, I.; Presnov, E.; Sprintsin, M.; Dag, A.
Evaluating water stress in irrigated olives: Correlation of soil water status, tree water status, and thermal
imagery. Irrig. Sci. 2009, 27, 367–376. [CrossRef]
79. Santesteban, L.G.; Di Gennaro, S.F.; Herrero-Langreo, A.; Miranda, C.; Royo, J.B.; Matese, A. High-resolution
UAV-based thermal imaging to estimate the instantaneous and seasonal variability of plant water status
within a vineyard. Agric. Water Manag. 2017, 183, 49–59. [CrossRef]
80. Espinoza, C.Z.; Khot, L.R.; Sankaran, S.; Jacoby, P.W. High Resolution Multispectral and Thermal Remote
Sensing-Based Water Stress Assessment in Subsurface Irrigated Grapevines. Remote Sens. 2017, 9, 961.
[CrossRef]
81. Baluja, J.; Diago, M.P.; Balda, P.; Zorer, R.; Meggio, F.; Morales, F.; Tardaguila, J. Assessment of vineyard
water status variability by thermal and multispectral imagery using an unmanned aerial vehicle (UAV). Irrig.
Sci. 2012, 30, 511–522. [CrossRef]
82. Li, H.; Luo, Y.; Zhao, C.; Yang, G. Estimating crop coefficients of winter wheat based on canopy spectral
vegetation indices. Nongye Gongcheng Xuebao/Trans. Chin. Soc. Agric. Eng. 2013, 29, 118–127. [CrossRef]
83. Zulini, L.; Rubinigg, M.; Zorer, R.; Bertamini, M. Effects of drought stress on chlorophyll fluorescence and
photosynthetic pigments in grapevine leaves (Vitis vinifera cv. ‘White Riesling’). Acta Horticulturae 2007, 754,
289–294. [CrossRef]
84. Zarco-Tejada, P.J.; Gonzalez-Dugo, V.; Berni, J.A.J. Fluorescence, temperature and narrow-band indices
acquired from a UAV platform for water stress detection using a micro-hyperspectral imager and a thermal
camera. Remote Sens. Environ. 2012, 117, 322–337. [CrossRef]
85. Gago, J.; Douthe, C.; Coopman, R.E.; Gallego, P.P.; Ribas-Carbo, M.; Flexas, J.; Escalona, J.; Medrano, H. UAVs
challenge to assess water stress for sustainable agriculture. Agric. Water Manag. 2015, 153, 9–19. [CrossRef]
86. Zarco-Tejada, P.J.; Guillen-Climent, M.L.; Hernandez-Clemente, R.; Catalina, A.; Gonzalez, M.R.; Martin, P.
Estimating leaf carotenoid content in vineyards using high resolution hyperspectral imagery acquired from
an unmanned aerial vehicle (UAV). Agric. For. Meteorol. 2013, 171, 281–294. [CrossRef]
87. Sagan, V.; Maimaitiyiming, M.; Fishman, J. Effects of Ambient Ozone on Soybean Biophysical Variables and
Mineral Nutrient Accumulation. Remote Sens. 2018, 10, 562. [CrossRef]
88. Urban, J.; Ingwers, M.W.; McGuire, M.A.; Teskey, R.O. Increase in leaf temperature opens stomata and
decouples net photosynthesis from stomatal conductance in Pinus taeda and Populus deltoides x nigra. J. Exp.
Bot. 2017, 68, 1757–1767. [CrossRef] [PubMed]
89. Zhang, J.; Cheng, D.; Li, Y.; Chen, H. Effect of Light and Water Stress on Photochemical Efficiency and
Pigment Composition of Sabina vulgaris Seedlings. Chin. Bull. Bot. 2017, 53, 278–289. [CrossRef]

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).