Land-productivity dynamics towards integrated assessment of land degradation at global scales.
Date: 2016-08-03
Institutions: Joint Research Centre
Subjects: biomass, data processing, degradation of the environment, desertification, environmental monitoring, land productivity, land use, remote sensing, research report

Land-Productivity Dynamics  
Towards integrated assessment of land 
degradation at global scales  

Ivits E., Cherlet M. 

2013 

EUR 26052 EN 

Reference time-seriesof vegetation indexReference time-seriesof vegetation indexReference time-seriesSavitzky Golay filter:savgol (N left, N right, degree)Nleft: 25Nright: 25degree of polynomial: 4Missing data handling: GAP > 2 decadesYES-> flagNO -> cubic spine interpolationIterative linear interpolationto days199920002001MI1 Std.1 Std.MINMINJ.D.SLE199920002001MI1 Std.1 Std.MINMINJ.D.SLEMI1 Std.1 Std.MINMINJ.D.SLEReference time-seriesof vegetation indexReference time-seriesof vegetation index3)365()365()365(200120001999SLESLESLEL 
 
 
 
 
 
This publication is a Technical report by the Joint Research Centre (JRC), the European Commission’s science and knowledge 
service. It aims to provide evidence-based scientific support to the European policy-making process. The scientific output 
expressed does not imply a policy position of the European Commission. Neither the European Commission nor any person acting 
on behalf of the Commission is responsible for the use which might be made of this publication. 

Contact information  
Name: Michael Cherlet 
Address: Joint Research Centre, Via Enrico Fermi 2749, 21027 Ispra (VA), Italy 
E-mail: Michael.cherlet@JRC.ec.europa.eu 
Tel.: +39 0332 789982 

JRC Science Hub 
https://ec.europa.eu/jrc 

JRC80541 

EUR 26052 EN 

PDF 

ISBN 978-92-79-32354-6 

ISSN 1831-9424 

doi:10.2788/59315 

Luxembourg: Publications Office of the European Union, 2016  

© European Union, 2016 

Reproduction is authorised provided the source is acknowledged. 

How to cite: Eva Ivits, Michael Cherlet; Land-Productivity Dynamics Towards integrated assessment of land degradation at global 
scales  title; EUR 26052; doi:10.2788/59315 

All images © European Union 2013 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1  Table of Contents 

1.  Background and Introduction ............................................................................................. 5 

2  Land-productivity dynamics map: Methodology................................................................. 7 

2.2  Derivation of phenological and productivity indices .................................................... 8 

2.3  Long-Term Change map of Land-Productivity .......................................................... 16 

2.3.1  The Steadiness Index: long term ecosystem change ............................................ 16 
2.3.2  Combining Steadiness Index with baseline levels for Standing Biomass ........... 24 
2.3.3  Standing Biomass state change ............................................................................ 28 
2.3.4  The Land-Productivity Long Term Change map ................................................. 30 

2.4  Current Status map of Land-Productivity. .................................................................. 32 

2.4.1  Global Phenological Types or Ecosystem Functional Types .............................. 32 

2.4.2 

Implementing the Local Net Scaling method ...................................................... 40 
2.5  Combined assessment of land-productivity ................................................................ 44 

3  Conclusion and Outlook .................................................................................................... 47 

Page 3 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Page 4 of 52 

 
 
 
 
 
1.  Background and Introduction 

The JRC Institute for Environment and Sustainability (IES) develops the new World 
Atlas  of  Desertification  (WAD)  in  support  of  European  Policies  addressing  regional  and 
global  Environmental  Issues,  Development  Cooperation  and  Integration  and  thus  working 
also  towards  successful  implementation  of  the  UN  Convention  to  Combat  Desertification 
(UNCCD). WAD contributes to  the Commission input in  the  UNCCD  process  and directly 
answers to the interest of the UNCCD stakeholders by providing a foundation to address the 
global  challenges  related  to  land  degradation  and  desertification.  WAD  is  considered  to 
become an international baseline to improve the implementation of UNCCD impact indicator 
reporting and to help design solutions for the implementation of the post Rio +20 process on 
obtaining  new  global  Sustainable  Development  Goals  (SDGs),  in  particular  the  target  on 
‘striving towards a land degradation neutral world’. 

To  address  these  complex  global  challenges,  a  monitoring  and  assessment  system 
offering up-to-date information on the status and trends of land degradation and their causes 
and effects is under development as routes for possible solutions need to be offered. A useful 
monitoring  and  assessment  system  will  supply  indicators  that  account  for  the  climate 
dependence of ecosystem functioning, is responsive to land cover and land use change while 
supplying  knowledge  of  the  temporal  and  spatial  patterns  of  ecosystem  dynamics  at  larger 
spatial  scales  (Ivits  et  al.,  2013a).  In  particular,  various  aspects  of  vegetation  productive 
dynamics  and  phenology,  reflecting  land  cover/use  transitions  that  can  lead  to  land 
degradation,  need  to  be  considered  in  a  spatio-temporal  context.  Because  of  the  large  areal 
coverage  and  continuous  temporal  sampling,  remotely  sensed  data  provides  a  synoptic 
picture  of  vegetation  dynamics  in  space  and  time  and  thus  have  a  great  potential  for 
monitoring  vegetation  and  ecosystem  change  from  regional  to  global  scales  (Myneni  et  al., 
1997).  

Building  on  numerous  studies  that  use  the  Normalized  Difference  Vegetation  Index 
(NDVI)  as  base  layer,  we  expand  this  set  of  variables  by  calculating  phenological  metrics 
from  time series of the  vegetation  index. By  de-convolution  of the original  time series into 
phenological  metrics  it  is  expected  to  yield  additional  information  on  various  aspects  of 
vegetation/land  cover  functional  composition  in  relation  to  dynamics  of  ecosystem 
functioning and land use (Ivits et al., 2013a). Using these vegetation phenological metrics can 
provide a quantitative basis to monitor such information on ecosystem dynamics and change. 
This is one of the central features in global change research as it is envisaged to provide users 
with  an  independent  measure  on  how  ecosystems  respond  to  external  impacts,  be  it  human 
induced or climate change (Linderholm, 2006; Parmesan, 2006; White et al., 2009, Fensholt, 
2012; Ivits et al., 2012b). 

The resulting remote sensing derived spatial layers are then combined with ancillary 
bio-physical and socio-economic information in order to flag areas that actually are affected 
by  land  degradation.  This  will  include  attributions  to  different  levels  of  intensity  and 
probability  of  major  causes,  which 
land 
in 
degradation/desertification issues as summarised by Sommer et. al., 2011 and shown in table 
1 below. 

include  major 

future  will 

the 

Page 5 of 52 

 
 
 
 
 
 
A.  Overuse of agricultural land, intensiﬁcation, inappropriate agricultural 

practices/non-SLM, increased soil erosion 

B.  Increase in intensive irrigation, overuse of water resources, salinization 

C.  Grazing mismanagement, overgrazing and decreasing NPP in rangelands, soil 

degradation, sand encroachment 

D.  Deforestation 

E.  Increased aridity or drought 

F.  Socio-economic issues, changes in population distribution 
and density, rural migration/land abandonment, urban 
sprawl 

G.  Uncontrolled expansion of mineral mining and industrial 
activities, extensive air and water pollution by waste 
materials, soil loss by contamination 

Table 1: Major desertiﬁcation and land degradation issues 

In the following chapters this report outlines the methodological details and the actual 
status of global and regional implementation of this new integrated assessment framework for 
land degradation at global scales. 

Page 6 of 52 

 
 
 
 
2  Land-productivity dynamics map: Methodology 

Conforming  to  the  concept  of  interacting  human-environment  systems,  it  is  human  activity 
and  demands  that  create  ‘ecosystem  services’.  Anthropogenic  impacts,  and  changes  over 
time, define the trends and the current potential for any ecosystem to supply these services. 
These demand-driven services tap into available ecosystem structures, which represent a vital 
natural resource, and affect the ecosystem functioning. The dynamics of the Earth’s covering 
biomass,  or  standing  biomass,  is  a  good  expression  of  the  general  level  of  the  potential  to 
supply, or keep on supplying, ecosystem services. The dynamics of the standing biomass in 
this study is evaluated by the change using time-series of long term observations. Assessing 
vegetative cover dynamics approximates a measure for general productivity levels of the land 
or  human-environment  system.  Land-productivity  reflects  climatic  constraints,  the  overall 
quality  of  the  land,  the  efficiency  of  using  the  land  and  other  resources  and  it  indirectly 
indicates the level at which these resources are appropriated for human use, i.e. land used for 
intensive/extensive  agriculture,  grazing,  forestry  or  urban  areas.  Hence,  land-productivity 
dynamics might indicate levels of sustained land-quality and is therefore used as first step in 
the land degradation assessment. Analysis of long-term changes and current efficiency levels 
of vegetative or standing biomass are combined into land-productivity dynamics according to 
the  scheme  presented  in  figure  1.  Efficiency  here  is  the  extent  to  which  the  biomass 
productivity  of  an  area  resembles  the  maximum  of  its  surrounding.  This  document  follows 
this flow.  

Figure 1: General scheme for the calculation of the land-productivity dynamics 

Page 7 of 52 

 
 
 
2.2  Derivation of phenological and productivity indices 

The derivation of phenological parameters from time-series of remote sensing images 
has  been  given  much  attention  (e.g.  Reed  et  al.,  1994;  Zhang  et  al.,  2003;  Jönsson  and 
Eklundh, 2004; Bradley et al., 2007). Most of these methods concentrate on the derivation of 
the  start,  the  end  and  the  maximum  date  of  the  vegetation  growing  season  together  with 
productivity  measures  approximating  Net  Primary  Productivity  and  growing  season 
productivity. However, for the characterization of ecosystem dynamics a more complete set 
of  indicators  is  needed.  Furthermore,  in  view  of  applying  the  method  on  the  continental  to 
global scale, a priori setting thresholds, whether or not adapted to local processes, should be 
avoided.  For  this  reason  the  “Phenolo”  software  was  developed  in-house  at  the  EC  Joint 
Research Centre.  

Phenolo envisages to derive a more complete range of variables to describe ecosystem 
dynamics at global scales. Phenolo therefore calculates all parameters and thresholds needed 
using the time series information for each single pixel. Hence, the software adapts to local as 
well as temporal variations and produces the most optimized set of variables for each pixel.  
The  following  gives  a  short  description  of  the  approach  while  Figure  2  presents  the  most 
important  steps  in  a  flowchart.  For  a  more  complete  description  of  the  method  and  the 
variables see Ivits et al., 2013a. 

The following Phenolo processing steps are explained below:   

Data preparation: 

(1) Data ingestion:  
(2) Data quality:  
(3) Pre-processing:  

rescaling of original data values 
gap filling and outlier adjustment 
temporal resampling and smoothing 

Processing:  

(4) time lag for moving averages 
(5) Intersections:  

phenological metrics 
productivity variables  

(1) The Phenolo software ingests time series data of varying length and from various 
satellite sensors (NOAA, SPOT, MERIS, etc.) and different products (NDVI, fAPAR, etc.). 
In order to render the calculated Phenolo variables comparable between the different satellite 
sensors,  products  and  temporal  ranges,  the  input  data  range  is  rescaled  and  values  are 
reported  according  the  interval  0%-100%  ([V0,  V1]).  In  case  in  the  input  dataset  the  0  is 
projected  to  the  value  0  the  rescaling  is  straightforward  and  the  output  range  is  defined  as 
V0=0 and V1=100 (i.e. [0,100]) in case of integer input values or V0=0 and V1=1 (i.e. [0,1]) in 
case of floating point values. In case in the input dataset 0 is not projected to the value 0 the 
data is rescaled accordingly where: 

0 = V0 – value in the dataset which corresponds to 0 and  
1 = (V1 –Vmin)/(Vmax-Vmin).  
For  instance,  in  case  of  an  input  dataset  with  a  physical  range  of  min=25  and 
max=275, V0 becomes 
0 = V0-25 thus V0 = 25 whereas 
1 = (V1-25)/(275-25) thus V1 = 275. 

Page 8 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
(2)  Some  time  series  have  missing  values  and  these  gaps  need  to  be  substituted  by 
acceptable values to obtain a continuous series to derive variables from. The GAP parameter, 
for which the Phenolo default is 2 decades, can be set according to the temporal resolution of 
the input time series. If the number of consecutively missing values in the time-series is less 
than the GAP parameter  these values  are replaced by  values  created through  a cubic spline 
interpolation.  Pixels  with  missing  time  steps  in  the  input  time-series  larger  than  the  GAP 
parameter  are  flagged  to  indicate  that  the  calculated  phenological  variables  are  uncertain. 
Negative outliers in the time series, that might be indicative for e.g. cloud cover, are removed 
by substituting each sufficiently pronounced downward spike in the time-series data with the 
cubic  interpolation  through  its  four  neighbors,  if  those  neighbors  are  contiguous;  if  they 
aren't, the value remains untouched. The number of outliers to be removed can be set by the 
user as well as the threshold for the outlier identification and subsequent removal. This latter 
value is given in percentage in proportion to the input NDVI range (rescaled between 0 and 
100, see before) and is calculated according to the following formula: 

(out-vb)*(out-va)<t2,                                                                                             (Eq. 1) 

where out is the value of the decade that is considered as outlier,  vb is the value before the 
outlier  data value,  va is  the value after the outlier  data value  and t  is  the value specified  in 
percentage of the input range. Thus, in order to define a decade (or any other time-step) as an 
outlier, the t parameter to the power of two has  to  be  larger  than the product  of the  outlier 
value  minus  the  time-series  value  before  and  the  outlier  value  minus  the  time-series  value 
after. 

(3) In order to generate results comparable between data sources with different time 
aggregation windows, the  next  pre-processing step involves interpolation  of the input time-
series temporal resolution to daily values. Another step smoothes the series through a linear 
interpolation.  Alternatively  a  spline  interpolation  can  be  chosen  for  the  calculation  of  the 
daily values. For both methods the number of iterations can be set by the user. Subsequently, 
an  iterative  Savitzky-Golay  filter  with  4th  polynomial  degrees  and  a  length  of  50  days  is 
applied  to  the  time  series  in  order  to  identify  and  remove  short  peaks  and  drop-offs  due  to 
noise as caused for instance by clouds. Both the window size and the degree of polynomials 
are parameters that can be adjusted according to best fit or need. These pre-processing steps 
(1-3)  result  in  the  reference  time  series  on  which  the  phenological  variables  are  to  be 
computed.  

(4)  The  methodology  for  calculating  the  phenological  and  productivity  variables  is 
based  on  an  original  method  by  Reed  et  al.  (1994)  and  uses  intersection  points  of  the 
reference time series with a backward and a forward shifted moving average smoothed time 
series to define two reference dates. For each pixel the forward and backward lag, i.e. the size 
of the moving average window, is defined by the length of the non-growing season. Reed et 
al. (1994) determined a pre-defined lag time based on their a-priori knowledge on the average 
phenology of the study area. However, considering that all global continents have a diversity 
of climatic regions, a large variations of ecosystems, miscellaneous land covers and a range 
of land uses, such a-priori  defined lag  cannot  account  for such diversity. Equally,  a pre-set 
threshold  to  indicate  the  reference  date  for  seasonal  dynamics,  as  commonly  used  e.g.  in 
Timesat  software,  cannot  account  for  such  diversity  either.  For  global  application  Phenolo 
needed to address the spatial and temporal specifics of a single pixel. Therefor the solution to 
calculate  the  time  lag  for  the  moving  average  series,  a  crucial  step  for  finding  meaningful 

Page 9 of 52 

 
 
 
 
 
intersection  points,  is  strictly  data  driven  and  uses  for  each  individual  pixel  its  own  time 
series dynamics to determine the lag.  

For  each  individual  pixel,  the  average  length  of  the  non-growing  season  that  will 
define the size of the moving average window, with which the original time series is lagged 
forward  or  backward  to  find  intersection  points,  is  calculated  using  the  yearly  estimated 
growing season length (SLE). This is derived from the pixel’s reference time-series by using 
two alternative methods (Figure 2).  

(a) 

(b) 

Between  two  subsequent  NDVI  signal  minima,  the  signal  above  the  line 
connecting  the  minimal  values  is  interpreted  as  a  histogram,  mean  and 
standard  deviation  of  which  are  determined.  The  obtained  mean  value  is 
effectively passing through the “barycentre” of the area, delimited by the 
NDVI signal and the baseline connecting the minima. By default, the SLE 
is  defined  to  be  two  times  the  standard  deviation  computed  from  the 
barycentre of the area  (expressed in  days).  Two  standard deviations were 
taken  as  the  default  value  for  SLE  because  that  describes  68.2%  of  the 
statistical population with normal distribution, which is considered a good 
approximation of vegetation productivity of the growing season.  
Alternatively, SLE is calculated as the side length of the quadrate, the area 
of which equals the area under the yearly vegetation index curve delimited 
by the two subsequent signal minima .  

The size of the moving average window, i.e. the lag, used to obtain the forward and 
backward shifted moving average filtered time series is then the yearly complement of SLE 
and is calculated as: 

                                                                                   (eq. 2) 

where  L  is  the  lag  (in  days),  N  is  the  number  of  years  in  the  time-series,  and  365  is  the 
number of days in the year.  

In this way, the time series dynamics of each pixel is incorporated in the derivation of 
the phenological metrics in an objective and user independent way, allowing the algorithm to 
be  applied  under  different  climatic  regions,  land-use  and  ecosystems.  Running  the  moving 
average in the forward direction (from the beginning to the end of the time series) results in a 
curve  that  lags  behind  the  reference  time-series  whereas  running  the  moving  average  filter 
backward creates a forward lagging curve.  

(5) The Season Begin Day (SBD) and Season End Day (SED) are determined as the 
intersections  of  the  reference  time-series  and  the  forward  and  backward  lagged  moving 
average  curves,  respectively.  Note  that  the  data  points  SBD  and  SED  cannot  always  be 
determined  unambiguously.  Determination  fails  when  there  is  no  significant  seasonal 
variation, either due to extremely low vegetation density, or because the vegetation cover is 
not  subject  to  seasonal  variations.  However,  the  failure  to  determine  these  points  in  itself 
allows  deducing  certain  land  cover  characteristics.  The  intersections  of  the  reference  time-
series  and  moving  average  curves  might  result  in  SBD  values  in  the  previous  year  and/or 
SED values in the next year. For the proper determination of these points in a given year the 

Page 10 of 52 

NSLELNi1)2365( 
 
 
 
  
 
 
 
 
time-series  should  contain  appropriate  number  of  decades  in  the  previous  and  in  the  next 
year. Once the SBD and SED data points are determined, a number of additional parameters 
can be computed which may allow to trace characteristics of the vegetation cover and their 
respective changes over time. Additional parameters comprise e.g. (see Figure 3): 

  season length: SL = SED-SBD  
  Maximum day and value: MXD and MXV 
  Standing  Biomass  (SB):  Approximation  of  gross  primary  production,  defined  as  the 

integral area comprising a+b+c+d+e+f+g (Minimum-Minimum Integral, MI). . 

  Background  fraction:  Minimum-Permanent  Integral  (MPI),  the  area  comprising 
a+b+c. This integral is not only characterising the presence of a perennial vegetation 
component  but,  depending  on  the  vegetation  index  used,  may  also  contain  a 
component of soil substrate.  

  Seasonal  vegetation  cover:  Cyclic  Fraction  (CF),  area  g.  This  parameter  can  be 

directly related to purely seasonal growth. 

  Permanent vegetation cover: Permanent Fraction (PF), area d+e+f. This parameter can 
be  related  to  the  vegetation  cover  that  does  not  have  a  seasonal  cycle  within  the 
growing year. 

  Season exceeding residual integral: SER, area d+f. This parameter can be assumed to 

be related the amount of senescent vegetation outside the growing season. 

Figure  5  presents  the  SPOT  VEGETATION  NDVI  signal  over  rainfed  arable  land 
pixels for nine ecozones in Europe (Figure 4) and the forward and backward lagged moving 
average curves which determine the SBD and SED points. The moving average curves adjust 
to the typical seasonality profile of each ecozone due to the lag which is shifted according to 
the time series dynamics of the different bio-climatic regions. The SBD point is not biased by 
false  season  starts  as  e.g.  in  the  northern  F  ecozone  and  the  moving  averages  are  not 
influenced by the strong inter-annual variability of NDVI minimums as e.g. in the H, I and J 
ecozones.  Figure  6  presents  examples  for  phenological  and  productivity  variables  averaged 
over  the  time-series  calculated  from  SPOT  Vegetation  time-series  data  using  Phenolo.  The 
productivity  variables  are  scaled  in  proportion  of  the  total  biomass  (MImn)  whereas  the 
phenological variables are expressed in days. 

Page 11 of 52 

 
 
 
 
Figure  2:  A  schematic  representation  of  the  main  Phenolo  processing  steps  for  a  time  series 
composed of three years. 

Page 12 of 52 

Reference time-seriesof vegetation indexReference time-seriesof vegetation indexReference time-seriesSavitzky Golay filter:savgol (N left, N right, degree)Nleft: 25Nright: 25degree of polynomial: 4Missing data handling: GAP > 2 decadesYES-> flagNO -> cubic spine interpolationIterative linear interpolationto days199920002001MI1 Std.1 Std.MINMINJ.D.SLE199920002001MI1 Std.1 Std.MINMINJ.D.SLEMI1 Std.1 Std.MINMINJ.D.SLEReference time-seriesof vegetation indexReference time-seriesof vegetation index3)365()365()365(200120001999SLESLESLELReferencetime-seriesBackwardMoving AverageForward Moving Average05010015020025019992000200120022003200420052006200720082009NDVIForward MABackward MA199920002001LLLLLL050100150200250200320042005050100150200250200320042005Jan.Dec.Calculation of integrals, dates and vegetation index valuesVegetation phenology and productivity indicators 
 
 
Figure 3: Schematic representation of the main phenological and productivity variables calculated by 
Phenolo (reprint from Ivits et al., 2013a). 

Ecozones: 

E: Cold and wet  

F: Extremely cold 
and mesic 

G: Cold and mesic 

H: Cool temperate 
and dry 

I: Cool temperate 
and xeric 

J: Cool temperate 
and moist 

K: Warm temperate 
and mesic  

L: Warm temperate 
and xeric  

N: Hot and dry 

Figure 4: European subset of the Global Environmental Classification with the nine ecozones 
(colours) and 73 strata (shown by numbered capital letters). 

Page 13 of 52 

Productivity parameters:SB=a+b+c+d+e+f+g (Standing Biomass)Approximation of the Total Surface Biomass.MPI=a+b+c (Minimum Permanent Integral)Area characterising the perennial vegetation component and depending on the vegetation index used may also contain components of soil substrate.SER=d+f (Season Exceeding Residual Integral)The amount of senescent vegetation  outside of the growing season. May be more indicative to before/after harvest of agricultural areas.CF=g (Cyclic fraction)This parameter may be directly related to purely seasonal growth.PF=d+e+f (Permanent Fraction)The amount of vegetation that does not have a characteristic seasonal cycle within the growing season.Phenological parameters:SBD = Season Begin Day.SL = Season Length.MXV = Maximum Value of the vegetation index. MXD = Maximum Day.SBV = Season Begin Value of the vegetation index.050100150200250200320042005050100150200250200320042005SLJan.Dec.abcdefgSEVMAXreference time-seriesof vegetation indexSBDSBVMBDMBVMEVMLMXDMXVSEDMED0%100% 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
E 

H

K

F

I 

L

G

J 

N

Figure 5: A six years subset of SPOT Vegetation NDVI profiles (solid lines) and the moving average 
curves (dashed: backward MA, solid grey: forward MA) under rainfed arable land pixels in the nine 
ecozones (reprint from Ivits et al., 2013a). 

Page 14 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figure 5: Examples for phenological and productivity variables calculated from the SPOT 
VEGETATION sensor with Phenolo averaged over the time-series. The productivity variables are 
scaled in proportion of the Standing Biomass (MImn) whereas the phenological variables are 
expressed in days. 

Page 15 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
2.3  Long-Term Change map of Land-Productivity 

2.3.1  The Steadiness Index: long term ecosystem change  

Most of the studies addressing ecosystem change rely on linear regression modelling 
of a monotonous trend in natural systems. This is disputable as it may not account well for 
the fluctuation inherent to natural systems (Ivits et al., 2013c). The methodology developed 
here refers to the fact that natural systems hardly change linearly or would react to human or 
natural  influences  in  a  linear  manner.  Therefore,  finding  a  systematic  model  that  describes 
spatially  distributed  areas  according  to  their  different  system  stages  in  terms  of  their 
likeliness  of  undergoing  ecosystem  change  or  their  rather  fluctuating  into  a  relative 
equilibrium  is  a  difficult  challenge.  Moreover,  the  linear  regression  as  a  quantitative 
parametric  model  relies  on  four  principal  assumptions  that  must  be  met.  If  any  of  these 
assumptions  are  violated  then  confidence  intervals  and  significance  tests  of  the  linear 
regression model may be misrepresentative. 

For  the  WAD  there  was  also  the  need  for  baseline  information  that  represents 
ecosystem  dynamics  during  a  rather  short  and  recent  time  window.  Linear  regression 
methods need to comply with strict statistical criteria, the violation of which prevents using 
the  significance  tests.  Also  due  to  the  variance  in  the  temporal  behaviour  of  the  observed 
ecosystem  the  scarce  significance  resulting  from  a  linear  regression  limits  enormously  the 
spatial  extent  of  the  sample  to  work  with.  For  land  degradation  as  well  as  for  other  global 
ecological studies the need of keeping a spatial continuous sample set outweighs the rigorous 
but limited possibilities of the linear regression.  

Ecosystem  change  therefore  is  described  in  our  method  as  the  long-term  tendency 
natural  systems  seem  to  move  towards  and  is  expressed  by  the  Steadiness  index.  The 
“Steadiness”  index  combines  the  trend  tendency,  i.e.  keeping  the  qualitative  notation  of  a 
positive or negative slope derived from linear regression combined with the net change of the 
time  series.  Defining  the  Steadiness  index  classes  the  linear  trend  (tendency)  and  the  net 
change  values  are  used  as  qualitative  classifiers  for  the  phenological  fluctuation  dynamics. 
Relying  on  a  convergence  of  evidence  rather  than  on  significance  values,  these  parameters 
are combined into a classification scheme that  expresses the long term direction into which 
the ecosystem appears to move in the given time window. The resulting classes represent the 
levels  of  steadiness  of  change  dynamics  of  the  phenological  or  productivity  metric.  These 
classes characterize ecosystems ranging from strong negative to strong positive dynamics and 
also  assign  classes  where  the  system  can  be  assumed  to  fluctuate  within  the  boundary 
condition  of  a  more  steady  equilibrium.  The  argument  is  that  a  not-significant  t-test  (or 
alternatively  the  non-parametric  Mann-Kendall  test)  of  the  linear  regression  should  not  be 
taken indicative of a non-changing ecosystem  state and that a simple measure like the here 
proposed  Steadiness  index  might  be  better  adapted  to  reflect  fluctuating  phenological 
dynamics of ecosystems. 

The Steadiness index is calculated by combining the following two simple measures:  

Page 16 of 52 

 
 
 
 
 
 
 
 
(1) The tendency of the change, expressed in the slope of the trend. Thus the tendency 
is calculated by fitting a linear trend over the time-series using the formula:  

                                                                                                         (eq. 3)  

  is  the  intercept, 

where 
  is  the  slope  of  the  fitted  line,  X  is  time  and  Y  is  the 
phenological  or  productivity  time  series.  The  slope  of  the  linear  trend  expresses  the 
dominant  tendency,  positive  or  negative,  towards  which  the  system  moves.  In 
assessing  this  tendency  no  hypothesis  test  is  performed  but  the  raw  values  are  used 
independently of significance to be integrated in a qualitative scheme. 

(2) The net change, positive or negative, of the phenological metrics over the selected 
time  period  supplies  a  second  evidence  of  ecosystem  dynamics  additionally  to  the 
slope  derived  from  the  linear  regression.  The  net  change  is  calculated  by  the  Multi 
Temporal  Image  Differencing  (MTID)  method  (Guo  et  al.  2008).  MTID  for  e.g  a 
Season Length time-series over the period 1982-2005 is calculated as: 

,                                                                                        

(eq.4) 

equals the Digital Number of the Season Length value in the corresponding 
where 
year.  Although  these  two  measures  develop  alike  over  large  areas,  due  to  the  non-
linear, fluctuating nature of Season Length values of ecosystems, a pixel with positive 
slope might express negative net change, and vice versa as shown in Figure 7.  

The  combination  of  the  tendency  and  change  of  the  time-series  provides  a 
convergence of evidence for monotonous and dominant,  upwards or downwards, change of 
the  system  and  results  in  the  classes  of  the  Steadiness  index.  There  are  four  potential 
combinations of the negative or positive trend and of the negative or positive change (Table 1 
and Figure 8) that represent the levels of equilibrium or possibly changing equilibrium of the 
phenological metric.  

The  Steadiness1  class  indicates  strong  negative  dynamics  during  the  observed  time 
series  with  monotonous  downwards  trends  and  negative  net  change,  indicating  that  the 
ecosystem is been developing towards a changing equilibrium.  

The Steadiness2 class shows areas where the trend  is negative but the net change is 
positive.  Here  the  system  fluctuates  but  does  not  show  a  clear  tendency  towards  a  new 
equilibrium characterised by e.g. systematically  shorter season length or  systematically less 
Standing  Biomass.  Hence,  these  areas  are  more  likely  to  remain  within  the  current 
equilibrium and have maintained a more steady state.  

The Steadiness3 class shows areas where the time series trends are positive but the net 
changes are negative. Here the system fluctuates in a different way but is likely to remain in 
current equilibrium and positive trends will not necessarily result in long term improvement 
or longer season of the ecosystems.  

Page 17 of 52 

XY1001200419822005)(iiDDMTIDiD 
 
 
 
 
 
 
 
 
The Steadiness4 class shows strong positive dynamics during the observed time series 
with monotonous upwards trends and net positive change,  indicating that the observed time 
series is probably changing equilibrium.  

These  four  classes  indicate  an  apparent  direction  into  which  the  phenological  or 
productivity dynamics of each pixel is evolving over time, upward or downward, confirmed 
or enhanced by the net change of the metric expressing the prevailing dynamic fluctuation of 
the  system.  The  method,  being  non-parametric,  does  not  have  to  comply  with  the 
assumptions  of  linear  regression  regarding  independence  of  the  errors,  stationarity  and 
normality  and  can  be  applied  on  any  type  of  data  or  time  window.  There  is  no  reliance  on 
statistical  significance  as  a  criterion  of  trend  relevance  in  order  to  avoid  the  necessity  of 
setting thresholds. Avoiding significance tests also enables the method to be applied on short 
time  series  where  the  calculation  of  statistical  significance  would  not  provide  meaningful 
information due to the limited number of observations. Instead, the Steadiness index is based 
on  a  convergence  of  evidences  that  the  ecosystem  changes  dynamics.  Furthermore,  this 
qualitative approach avoids the need to introduce possibly subjective or solely locally valid 
thresholds  for  interpretation  of  strengths  of  trends  and  is  not  spatially  restrictive  based  on 
subjective thresholds.  

Steadiness classes: 

STEADINESS1: negative slope and negative change.  

Represents pixels under strong and negative ecosystems dynamics, with a 
probability of changing equilibrium. 

STEADINESS2: negative slope and positive change.  

Represents pixels under moderate negative ecosystems dynamics but likely to 
remain in current equilibrium.  

STEADINESS3: positive slope and negative change.  

Represents pixels under moderate positive ecosystems dynamics but likely to 
remain in current equilibrium  

STEADINESS4: positive slope and positive change.  

Represents pixels under strong and positive ecosystems dynamics with a 
probability of changing equilibrium. 

Table 1: Summary of the four Steadiness index classes 

Page 18 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Steadiness1: Examples of negative slope (linear trend) and negative 
change (MTID) of the SL time-series for the three zones 

Steadiness2: Examples of negative slope (linear trend) and positive 
change (MTID) of the SL time-series for the three zones  

Steadiness3: Examples of positive slope (linear trend) and negative 
change (MTID) of the SL time-series for the three zones 

Steadiness4: Examples of positive slope (linear trend) and positive 
change (MTID) of the SL time-series for the three zones 

Figure 7: Slope (upper left) and change (MTID, upper right) of the linear trend of the Season Length 
(SL)  time-  series  for  the  years  1982-2005.  Graphs:  The  SL  time-series  averaged  within  the  zones 
North,  Central  and  South  as  in  Figure  3.  (Y  axis:  Season  Length  in  days;  X  axis:  calendar  year). 
Boxes indicate the calculated linear trend equation respectively the MTID value of the presented SL 
time series. Reprint from Ivits et al., 2013c. 

It  must  be  noted  that  Steadiness,  as  much  as  linear  regression,  cannot  model  those 
ecosystem  changes  that  follow  a  non-linear  pattern.  However,  for  most  ecosystem  change 
studies  the  main  interest  is  not  the  form  a  non-linear  change  takes  but  rather  the  general 
direction  the  system  evolves  to.  Therefore,  relying  on  convergence  of  evidence  from  the 
slope of the linear  regression and from  the change  (MTID) indicator of the time-series, the 
steadiness replaces the significance test and assesses the general nature of the change without 
excluding pixels from the further analysis  

Furthermore,  we  note  that  the  last  value  of  the  time  series,  which  is  used  as  the 
reference in the MTID index, might have a strong effect on both the calculated MTID value 

Page 19 of 52 

 
 
 
 
 
 
 
 
and on the slope value in case the least  year is  an outlier. This  is  not  an entirely unwanted 
effect for two reasons.  

(1)  In  our  study  we  search  for  convergence  of  evidence  for  similarity  in  the 
change of ecosystem dynamicss therefore it is desirable to use change measures that 
are affected in the same way. Using e.g. eq. 2 in Guo et al. (2008) would introduce a 
change  indicator  of  another  nature  which  would  disable  finding  convergence  of 
evidence  because  of  showing  ecosystem  changes  in  a  different  manner.  The  MTID 
measure of equation 2 in Guo’s paper takes the absolute value of each paired change 
values i.e:  

                                                                                          (Eq. 5) 

and  therefore  measures  the  change  intensity  indicator,  which  is  an  entirely  different 
issue and thus cannot be used for finding convergence of evidence. One could change 
the equation and not take the absolute value of the paired changes in the form of: 

                                                                                          (Eq. 6) 

Please note however that in this case the equation mathematically equals the last year 
minus  the  first  year  of  the  time-series.  Therefore  that  measure  would  be  even  more 
affected by the  first  and/or the last  year than  eq.  1 what  we have used in  this study 
but, unlike eq. 1, does not take the intermediate years into account.  

(2)  At  such  early  stage  of  a  time-series  analysis  it  is  not  possible  to  state 
weather the first or the last year is a real outlier or rather a measurement error e.g. a 
sensor artefact. Change measures like e.g. Eq. 5 and 6 have a strong smoothing effect 
both on the magnitude of the changes as well as on the effect of the first and/or last 
values in  the time-series.  If extreme values in  the first  or last  years are not  outliers, 
these values form an important component of the time series, as an extremely wet or 
extremely dry year or significant land use change might have happened just in one of 
those  years.  Therefore,  smoothing  out  the  effect  of  these  values  at  this  stage  of  the 
analysis could be a mistake that might lead to severe information loss. Therefore, in 
this  study  eq.  4  is  used  for  assessing  the  change  while  being  well  aware  of  the 
possible effect of the last year on the value and we argue that this affect can only be 
accounted for with extra information on land use change and climatic effects.  

Page 20 of 52 

200419821iiiDDMTID200419821iiiDDMTID 
 
 
 
 
 
 
 
 
 
Steadiness Classes for 
Standing Biomass 

Steadiness Classes for 
Standing Biomass 

Figure 8: A global and a European example for the Standing Biomass Steadiness classes derived from the 
NOAA GIMMS3g dataset from 1982-2010. 

Page 21 of 52 

 
 
 
Calculation of Land-Productivity Long Term Change Map 

The  Steadiness  Index  provides  a  good  measure  for  the  general  track  the  productive 
capacity  of  an  ecosystem  is  moving  on.  This  knowledge  is  very  useful  for  performance  or 
change analysis as it gives a measure on the progression of a general process. The term land 
‘degradation’ certainly implies that a process needs to be assessed. But in order to estimate 
the  importance  of  an  on-going  process  in  terms  of  ‘degradation’  one  needs  to  know  the 
baseline  from  where  the  process  started  and  also  how  strong  this  process  is  or  was.  In  the 
framework  of  land  degradation  and  desertification  assessment,  the  Steadiness  Index  alone 
does not provide all the information that is needed to qualify the observed process. One needs 
to  know  as  well  whether  an  ecosystem  changed  its  state,  see  Figure  9,  and  a  negative  or 
positive Steadiness Index does not necessarily provide that information.  

Therefore, further analysis will need information added that (1) quantifies the original 
productivity  level  where  the  system  comes  from  (2)  quantifies  the  strength  of  the  negative 
and positive change, determining if e.g. the Standing Biomass level of the ecosystem changed 
class  during  the  period  of  interest.  This  assesses  where  the  ecosystem  arrived  to  after  e.g. 
negative  change  dynamics.  Even  though  an  ecosystem  expressed  long  term  negative 
steadiness and it was shown that this negative change was strong in certain high productivity 
systems this long term negative change might not be enough to for the system to change its 
equilibrium  and  arrive  to  a  lower  productivity  state  (Figure  10).  Similarly,  if  the  system 
shows positive dynamics but that dynamics was not strong enough the system still might be 
in low productivity conditions. 

Figure  9:  Ecosystem  state  change  when  resilience  thresholds  are  passes  (After  S.  Prince, 
2012 Personal Communication).  

Page 22 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
Figure 10: Schematic representation of Standing Biomass change classes of ecosystems . The y-axis 
represents  the  levels  of  Standing  Biomass  whereas  the  x  axis  represents  time.  St1-4  stands  for  the 
Steadiness classes. WK=weak; MOD=moderate; STR=strong. 

Page 23 of 52 

 
 
 
 
 
 
2.3.2  Combining Steadiness Index with baseline levels for Standing Biomass 

Figure  10  showed  schematically  the  importance  to  establish  the  baseline  start  to 
determine the relative importance of the ecosystem change dynamics that have been captured 
by the Steadiness. This next analysis was performed in order to arrange the levels of Standing 
Biomass at the beginning of the time series into three categories: low, medium and high. In 
order to quantify the levels of Standing Biomass at the beginning of the time-series, Standing 
Biomass  values  for  the  first  three  years  of  the  series  (i.e.  1982-1984)  were  temporally 
averaged. Averaging fewer  years  would bias the estimation  of the initial  Standing  Biomass 
classes  because  of  the  influence  of  extreme  values  due  to  e.g.  drought  of  very  wet  years. 
Averaging more  years  for the quantification of the  Standing Biomass  levels  would bias the 
assessment  towards  the  series’  mean  value.  In  order  to  arrange  the  scale  values  into 
categories ISODATA (Iterative Self Organising Data Analysis technique) classification was 
run on the three years temporal average. The choice of unsupervised classification was opted 
for  because  of  lack  of  proper  data  or  literature  source  relating  satellite  measured  Standing 
Biomass estimates to field measured Standing Biomass or NPP values at the global scale.  

the  Worlds 

The ISODATA algorithm arranges the input data population (i.e. the number of land 
surface pixels) into homogeneous clusters where the clusters have close to equal number of 
members.  Therefore,  when  running  the  classification  10  classes  were  opted  for  in  order  to 
produce categories such way that each of them contains ca. 10% of the data. Approximately 
40%  of 
land  resources  are  drylands  (Drylands  Development  Centre, 
www.undp.org)  which  can  be  assumed  to  represent  40%  of  the  Globe’s  low  Standing 
Biomass areas. Therefore, the first four classes of the ISODATA run were assigned to be low 
Standing Biomass classes. The five consecutive classes were assigned to be medium category 
whereas the upper 10% of the values was selected to represent areas with high productivity 
levels. In case the input image were classified into three categories only, the input population 
would be homogeneously distributed among the three categories and differences between the 
Standing Biomass levels would not be properly addressed. Figure 11 shows the low, medium 
and  high  initial  Standing  Biomass  levels  classified  directly  into  three  categories  (Figure  11 
top) as well as the reclassification from 10 ISODATA classes (Figure 11, bottom).  

SB classes 

minimum 

maximum 

Class1 
Class2 
Class3 
Class4 
Class5 
Class6 
Class7 
Class8 
Class9 
Class10 

0 
402 
664 
840 
982 
1123 
1261 
1407 
1613 
1945 

402 
664 
840 
982 
1123 
1261 
1407 
1613 
1945 
2990.3 

mean 

239.6 
542.5 
756.7 
908.8 
1051.1 
1195.1 
1333.2 
1505.3 
1770.5 
2175.3 

Standard 
deviation 
86.4 
77.7 
51.1 
40.6 
40.5 
39.9 
41.8 
59.7 
94.7 
175.1 

Table 2: Class statistics for the initial Standing Biomass (SB) classes averaged for 1982-1984 

Page 24 of 52 

 
 
 
 
 
 
 
 
 
 
Initial Standing 
Biomass Levels 

Initial Standing 
Biomass Levels 

Figure 11: ISODATA classification of the initial Standing Biomass levels where the initial Standing 
Biomass levels represent the average values of the first three years (1982-1984). 

Page 25 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
To obtain a single representation that combines the Steadiness Index, indicating the general 
track that the ecosystem change is on-going, with the baseline levels of productive capacity 
expressed by three categories of Standing Biomass averaged during the first three years of the 
satellite time series, a lookup table is used as shown in table 3. The main colours used in table 
3 and the combined class annotations (e.g. St1, low) are used on the maps shown in figure 12.  

Steadiness classes 

Initial SB levels 

Steadiness 1 

Steadiness 2 

Steadiness 3 

Steadiness 4 

Low 

Medium 

High 

1      St1, low 

4           St2, low 

7         St3, low 

10     St4, low 

2      St1, medium  5           St2, medium  8         St3, medium  11     St4, medium 

3      St1, high 

6           St2, high 

9         St3, high  

12     St4, high 

Table 3: Look up table for the combination of the steadiness classes with the initial Standing Biomass 
(SB) levels. St1-4 stands for the 4 Steadiness classes. 

Page 26 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Steadiness Index 
combined with 
baseline levels of 
Standing Biomass 

Steadiness Index combined 
with baseline levels of 
Standing Biomass 

Figure  12:  A  global  and  a  European  mapping  of  the  combination  of  the  Steadiness  Index  and  the 
baseline levels of Standing Biomass. 

Page 27 of 52 

 
 
 
 
 
 
 
 
2.3.3  Standing Biomass state change 

Further  to  the  general  direction  in  which  the  ecosystem  productive  capacity  is 
moving,  expressed  by  the  Steadiness  Index,  and  the  initial  level  from  which  the  dynamics 
moved, for degradation assessment it is important to know whether the environment system 
also changed productivity state. This is important because when productivity thresholds are 
passed the environment system reaches another state. Natural or human induced pressures or 
opportunities  could  cause  this  state  change.  The  state  change  reflects  either  that  natural 
resilience thresholds have been passed or that other land uses/practices were introduced.  

In order to calculate the change classes for the Standing Biomass (SB), further to the 
above established classes at the start of the observation time period (Table 3 and Figure 12), 
the SB levels at the end of the time series were also computed. First, the last three years of 
the series (i.e. 2008-2010) SB values were temporally averaged. Second, the averaged values 
were submitted to an ISODATA classification in order to categorise the SB values at the end 
of the time series.  Using  a  classification  with  three  classes  at  the beginning  and  end of the 
time  series  would  only  allow  pinpointing  very  large  changes  especially  in  the  low  and 
medium classes due to the size of these classes. This means that the SB value would need to 
change considerably before changing from the medium to lower  productivity, or vice versa. 
Whereas, due to the relative small and more condensed value range of the high SB class,  a 
small  decline  might  already  indicate  a  significant  productivity  change.  Therefore,  the 
ISODATA  classification  with  10  classes  was  used  at  both  ends  of  the  time  frame.  The 
thresholds  of  the  10  ISODATA  classes  derived  from  the  initial  three  years  (Table  2)  were 
applied  similarly  over  the  last  three  years  average  SB  images  and  10  end  SB  classes  were 
assigned accordingly. This  way the comparison  of the  start  and  end SB  levels to  determine 
class change is based on the initial SB baseline classes. The initial and end 10 SB levels were 
compared and three change classes were assigned: (1) the NPP value remained in the initial 
class; (2) the NPP value changed 1 class; (3) the NPP value changed 2 or more classes (table 
4 and figure 14).  

Standing Biomass Class at END of time series  

Standing 
Biomass 
at 
START 
of time 
series 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 
2 
3 
4 
5 
6 
7 
8 
9 
10 

1 
2 
3 
3 
3 
3 
3 
3 
3 
3 

2 
1 
2 
3 
3 
3 
3 
3 
3 
3 

3 
2 
1 
2 
3 
3 
3 
3 
3 
3 

3 
3 
2 
1 
2 
3 
3 
3 
3 
3 

3 
3 
3 
2 
1 
2 
3 
3 
3 
3 

3 
3 
3 
3 
2 
1 
2 
3 
3 
3 

3 
3 
3 
3 
3 
2 
1 
2 
3 
3 

3 
3 
3 
3 
3 
3 
2 
1 
2 
3 

3 
3 
3 
3 
3 
3 
3 
2 
1 
2 

3 
3 
3 
3 
3 
3 
3 
3 
2 
1 

Table 4: Look up table for the classification of Standing Biomass change classes 

1: no change; 2: changed 1 class; 3: changed 2 or more classes 

Page 28 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
The maps in figure 13 illustrate the 10 classes of Standing Biomass levels at the start of the 
observation period with values averaged for the years 1982-1984 (top) , and at the end of the 
series averaged over 2008-2010 (bottom).   

Class level of 
Standing 
Biomass at 
START of time 
series (initial) 

Class level of 
Standing 
Biomass at 
END of time 
series (end) 

Figure 13: ISODATA classification of the Standing Biomass levels at the beginning and at the end of 
the  time  series  applied  over  the  temporal  average  of  three  years  (1982-1084  and  2008-2010, 
respectively). 

Page 29 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Class Change for the 
Standing Biomass variable 
over the period 1982-2012 

Figure 14: Class Change of Standing Biomass using the look up table shown in table 4. 

2.3.4  The Land-Productivity Long Term Change map 

The  ecosystem  productivity  Long  Term  Change  map  is  one  of  the  two  input  layers 
together with the ‘current status map’, that will be described further below, that are the basis 
for the final Land-Productivity Dynamics map. The Long Term Change map is obtained by 
combining  the  data  layers  derived  before:  ‘the  Steadiness  plus  initial  Standing  Biomass 
levels’ data with the ‘Class Changes of Standing Biomass’. The lookup table applied for this 
combination is shown in table 5 and resulting maps are illustrated in figure 15. . 

St1 
low 

St1 
med. 

St1 
high 

St2 
low 

St2 
med. 

St2 
high 

4 

10 

5 

11 

6 

12 

St3 
low 

7 

13 

St3 
med. 

St3 
high 

8 

14 

9 

15 

St4 
low 

10 

16 

St4 
med. 

St4 
high 

11 

19 

12 

22 

10 

11 

12 

13 

14 

15 

17 

20 

22 

10 

11 

12 

13 

14 

15 

18 

21 

22 

1 

1 

2 

3 

2 

4 

5 

6 

3 

7 

8 

9 

No 
change 

Changed 
1 class 

Changed 
2/more 
classes 

Table  5:  Look  up  table  for  the  combination  of  the  Steadiness  and  initial  NPP  levels  with  the  NPP 
change classes (see section 2.3.3). St1-4 stands for the Steadiness classes; med=medium. 

Page 30 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Land-Productivity Long Term Change  
(Standing Biomass period 1982-2010) 

Figure  15:  A  global  and  a  European  example  of  the  Long  Term  Change  Map  where  the  Standing 
Biomass  (SB)  Steadiness  Index  is  combined  with  the  baseline  SB  and  with  the  Class  Change  of SB 
applying the classification scheme shown in Table 5. 

Page 31 of 52 

 
 
 
 
 
2.4  Current Status map of Land-Productivity.  

Localized and contextual measurements of current levels of productivity efficiency are 
needed  as  additional  information  to  relate  long term  productivity dynamics  to  potential  on-
going  land  degradation  processes.  The  Local  Net  Scaling  (LNS)  approach  (Prince  2004, 
Wessels et al., 2007) answers this need where remote sensing estimated productivity of each 
pixel is expressed relative to the 90 percentile observed in all pixels falling within the same 
homogeneous environmental unit. Compilation of an adequate stratification of homogeneous 
biophysical units is crucial for the LNS method. The below section explains the compilation 
of  this  stratification  and  the  implementation  of  the  LNS  method  to  derive  the  land-
productivity  status  map  that  is  a  second  base  layer  needed  to  compose  the  final  land-
productivity dynamics map.  

2.4.1  Global Phenological Types or Ecosystem Functional Types 

The  LNS  method  requires  the  stratification  of  environmental  units  with  similar 
production potential, which is defined by spatial information on vegetation, soils, terrain and 
climate  (Wessels  et  al.,  2008).  Besides  this  biophysical  information,  ecosystems  may  be 
further characterized by  the physiognomy and functional  dynamics of the vegetation  cover. 
When characterising ecosystems,  as for  application of the  LNS  method, a compound  set  of 
functional attributes that describe vegetation dynamics should also be included.  

Phenological  types  are  defined  as  spatial  units  with  similar  patterns  of  seasonal 
phenology  and  productivity  dynamics  which  exhibit  similar  responses  to  changing  natural 
and human induced environmental conditions following the ideas of Paruelo et al. (2001) and 
Stow et al. (2000). Ivits et al. (2013a,d) described the compilation in detail and illustrated that 
such spatial units, named Ecosystem Functional Types (EFTs), reflect both climate and land 
use situations on the continental (Ivits et al, 2013a) and global level (Ivits et al., 2013d) and 
therefore  offer  a  meaningful,  transparent  and  objective  stratification  that  supports  applying 
the LNS method as one of the cornerstones for global land degradation studies. 

The  computed  Phenolo  parameters  were  summarized  in  their  temporal  mean  (1982-
2010)  and  consecutively  screened  against  multicollinearity  based  on  the  correlation  matrix. 
Variables  with  very  high  correlation  (>0.7)  were  removed  from  the  analysis.  Principal 
Component Analysis (PCA) was run on the correlation rather than on the covariance matrix 
in  order  to  standardize  the  input  Phenolo  variables  with  different  measurement  scales.  The 
Principal Components were rotated with the “varimax” technique and post-normalized with 
the  Kaiser’s  method.  The  rotation  of  the  PCA  eigenmatrix  was  chosen  in  order  to  clearly 
associate each PC axes with one Phenolo variable. This first, "screening PCA" served for (1) 
the selection of that set  of Principal Components that explained the highest amount of total 
variance  in  the  Phenolo  variables  measured  along  the  global  ecosystems  and  (2)  for  the 
selection  of  those  Phenolo  variables  that  demonstrated  the  highest  loadings  on  the  selected 
components. The loadings of the individual variables were normalized by multiplication with 
the square root of the eigenvalues in order to present the values as correlation with the PCA 

Page 32 of 52 

 
 
 
 
 
 
 
 
 
axes. With rotating the factor loadings the same amount of total explained variation could be 
reached  with  the  first  5  PCA  axes  but  the  explained  variation  was  more  homogeneously 
distributed  among  the  rotated  axes  (Table  6).  Furthermore,  with  the  rotation  each  PCA 
dimension  became  clearly  associated  with  one  variable  only  (Table  7).  MI  was  selected  to 
represent the first PC axis because the maximum value (MXV) might be a subject of outliers 
even after smoothing the vegetation index time series as explained in section 2.2. Based on 
highest  loadings  on  the  other  rotated  PC  axes  the  Cyclic  Fraction,  the  Maximum  Day,  the 
Season  Length  and  the  Permanent  Fraction  were  selected  to  enter  the  follow  up  analysis 
(table 7). Once the final number of principal  components was determined and the variables 
with the highest loadings were selected a final rotated PCA was run to clearly associated each 
pre-selected  Phenolo  variable  with  one  PC  axis  (Table  8).  The  spatial  pattern  of  the  final 
rotated  components,  i.e.  the  Empirical  Orthogonal  Functions  (EOFs),  was  calculated  by 
multiplying the loadings with the selected phenological variables as shown in figure 16.  

Table 6: PCA of the phenological variables and the global Phenological Types clusters. Statistics for 
the first five initial components and for the Varimax rotated components are shown. 

Table 7: Normalized loadings of the phenological variables on the first five original (left) and rotated 
(right) PCA axes. 

Rotated PCA components 
PC1 
0.091 
0.212 
-0.005 
0.947 
0.373 

PC2 
-0.043 
0.952 
-0.057 
0.212 
0.360 

CF 
MI 
MXD 
SL 
PF 
Table  8:  Normalized  loadings  of  the  selected  phenological  variables  on  the  first  five  rotated  PCA 
axes in the final PCA analysis. 

PC5 
0.188 
0.206 
-0.028 
0.218 
0.778 

PC4 
0.123 
-0.065 
0.992 
-0.002 
-0.052 

PC3 
0.969 
-0.048 
0.108 
0.096 
0.351 

Page 33 of 52 

 
 
 
  
 
 
 
 
 
A

B

C

D

E

Figure  16:  The  first  five  rotated  Empirical  Orthogonal  Functions  calculated  with  a  PCA  from  the 
phenological variables SB, Cyclic Fraction, MXD, SL and PF. 

ISODATA  cluster  analysis  was  run  on  the  rotated  EOFs  calculated  from  the  final 
PCA  model.  The  cluster  analysis  was  run  with  5  iterations,  min=10  and  max=500  classes 
requested, and the maximum class standard deviation set to 4 in order to allow homogeneous 
clusters. The clustering process was run on the EOFs rather than on the original phenological 
variables  because  through  the  PCA  with  correlation  matrices  the  resulting  eigenvectors  are 

Page 34 of 52 

 
 
 
 
 
 
 
normalised with zero mean and 1 SD. This delivers statistical features of the eigenvectors that 
are desirable for the calculation  of class means evenly distributed in  the data space and for 
the  iterative  clustering  of  the  pixels  using  minimum  distance  techniques.  The  resulting 
homogeneous ISODATA clusters represent major Ecological Functional Types (EFTs) of the 
global ecosystems (see figure 17).  

Figure 17: Top: Isodata classification of the first 5 rotated PCA axes. The classification resulted in 
100  homogeneous  EFTs.  Bottom:  Hierarchic  Isodata  classification  of  the  100  PTs  with  the  rotated 
EOF values averaged within the clusters. The classification resulted in 14 homogeneous EFTs.  

In order to ease the interpretation of the 100 EFTs the spatial units were submitted in 
another hierarchic ISODATA clustering process. Within each of the 100 EFTs the mean EOF 
values  of  the  five  selected  dimensions  were  calculated  and  written  in  a  new  multi-band 
image. This new image with the 5 averaged EOF values entered a new ISODATA cluster run. 
The setting of the clustering process was as before but the maximum numbers of classes were 
set to 100. This created a new, hierarchical grouping of the EFTs. The process of averaging 
the  EOF  values  within  the  new  clusters  were  repeated  as  long  as  the  number  of  resulting 
clusters converged and stabilized. All together two more clustering processes were necessary 
and the number of clusters stabilized at 14 (Figure 17, bottom). 

In  order  to  correlate  the  EFTs  map  of  global  ecosystems  with  existing  datasets 
Detrended Correspondence Analysis (DCA) was carried out with the Köppen Climate Zones 

Page 35 of 52 

 
 
 
 
 
 
(Peel at al. 2007) and the FAO land use system classes (FAO, 2011) . The 15 FAO classes 
were  represented  according  to  the  climate  zone  Tropical  (Tr),  Arid  (Ar),  Temperate  (Tm), 
and Cold (Cd) they are located in. The CA axes were detrended by 2nd order polynomials as 
suggested by Ter Braak and Prentice (1987) in order to avoid the arch effect (Hill and Gauch, 
1980)  caused  by  the  strong  environmental  gradients  and  for  better  interpretability  of  the 
ordination bi-plot. The DCA was carried out over the cross-tabulation of the EFTs with the 
Köppen classes and with the FAO classes, respectively. Results were presented in a bi-plot of 
the  first  two  dimensions  where  close  Euclidean  distances  between  the  spatial  units  of  the 
EFTs  and  the  Köppen  classes  respectively  the  FAO  land  use  classes  indicate  good 
correspondence with the phenology driven classification. Spatial units that lie at the centre of 
the  ordination  biplot  are  not  well  represented  by  the  dimensions  derived  from  the 
phenological data whereas spatial units at the edge of the biplot represent areas with very few 
observations. In order to demonstrate those EFTs that are not well represented, i.e. where the 
correspondence  to  the  Köppen  Zones  respectively  the  FAO  classes  was  low,  the  absolute 
values  of  the  sample  scores  were  added  over  the  first  two  dimensions.  The  sample  scores 
represent the distance of each sample from the centre of the biplot, therefore low cumulative 
sample scores over the first two dimensions indicate areas that are not well represented by the 
DCA biplot (the ones that are plotted in the center of the biplot). These values were mapped 
for  both  DCA  analyses.  Note,  this  is  not  a  general  goodness  of  fit  measure;  it  merely 
represents  those  areas  that  are  not  strongly  related  to  the  first  two  DCA  dimensions.  Other 
dimensions might represent these areas better, however given their lower explanatory power 
those associations are much less significant. 

The first four DCA axes explained 59.8% of the variation in the data, showing good 
correspondence between the Köppen Zones and the ISODATA classes (table 8). The first two 
axes explained 48.8% of the variance and as shown by the bi-plot most of the Isodata clusters 
could  be  associated  to  one  of  the  Köppen  Zones.  For  instance,  the  Dsd,  Dwd  and  Dfd 
(Cold/dry  summer/cold  winter,  Cold/dry  winter/cold  winter  and  Cold/no  dry  season/cold 
winter, respectively) zones were associated with the ISODATA clusters 40, 41, 47 with late 
maximum, high Cyclic Fraction, low Standing Biomass with short season and low PF. On the 
other  end  of  the  first  axis’  gradient  the  BWh  zone  (Arid/Desert/Hot)  and  the  Cwc  zone 
(Temperate/Dry  winter/Cold  summer)  was  associated  to  the  clusters  1,  5,  7  with  moderate 
SB, low PF and Cyclic Fraction, early maximum, and short seasons whereas the Köppen zone 
BWk  (Arid/Desert/Cold)  is  more  associated  to  the  clusters  30-32  with  lowest  SB  and  PF, 
short season, high CF and late maximum. On the far end of the second DCA axis the zones 
Af  (Tropical/Rainforest)  and  Am  (Tropical/Monsoon)  were  strongly  associated  with  the 
clusters 18, 55 and 58 with high SB, moderate to high season length and Permanent Fraction, 
low Cyclic Fraction and early maximum. Figure 18 presents for each EFTs the fit of first four 
axes of the DCA ordination with the Köppen Zones. Most ecosystems present a good (>50%) 
or  very  good  (>70%)  fit  with  the  ordination  axes  showing  a  good  fit  between  the  Köppen 
Zones  and  the  EFTs.  Notable  disagreements  were  observed  over  the  Great  Plains  of  the 
United States,  over continental Europe, the central  part of the Sahel  and  along the Yangtze 
river basin in China (figure 19). 

Detrended Correspondence Analysis (DCA) of the EFTs with the Köppen climate zones. 

Eigenvalues 
Cumulative % of explained variance 

Axis1 
0.647 
37.3 

Axis2 
0.201 
48.8 

Axis3 
0.145 
57.2 

Axis4 
0.060 
59.8 

Table 9: Detrended Correspondence Analysis (DCA) of the Isodata classes with the Köppen climate 
zones. 

Page 36 of 52 

 
 
 
Figure  18:  Detrended  Correspondence  Analysis  triplot  of  the  EFTs  clusters  classified  in  14 
categories  (icons  and  colors)  and  the  Köppen  climate  zones  (crosses).  The  zone  Cwc  (see  Table 
below)  was  excluded  from  the  analysis  due  to  low  number  of  samples.  Arrows  represent  the 
phenological  and  productivity  variables  which  are  passively  projected  onto  the  axes  thus  indicate 
correlation to the gradient. 

Köppen climatic zones and FAO Land Use Systems classes used in the present study 
LUS classes 
Forests, unmanaged 
Forests, unmanaged 
Grassland, 
unmanaged 

Köppen –Geiger 
Tropical rainforest 
Tropical Monsoon 

Code  Köppen -Geiger 
Af 
Am 

Code 
FRu 
FRm 

Tropical Savannah 

GRu 

Aw 

Cold, dry and warm summer 

Temperate, dry winter, cold summer 
Cold, dry and hot summer  

SHu 

SHm 

Grassland, managed  GRm 
Shrubland, 
unmanaged 
Shrubland, 
managed 
Rainfed agriculture 
AGr 
Irrigated agriculture  AGi 
Sparse  vegetation, 
unmanaged 
Sparse  vegetation, 
managed 
Wetlands 

WTL 

SPm 

SPu 

Arid, desert, hot 
Arid, desert, cold 
Arid, steppe, hot 

Arid, steppe, cold 

Temperate, Dry and hot summer  

Temperate, Dry and warm summer 
Temperate, no dry season, hot summer 
Temperate,  no  dry  season,  warm 
summer 
Temperate,  no  dry 
summer 
Temperate, dry winter, hot summer 
Temperate, dry winter, warm summer 

season,  cold 

Bwh 
Bwk 
Bsh 

Bsk 

Csa 

Csb 
Cfa 

Cfb 

Cfc 

Cwa 
Cwb 

Cold, dry and cold summer 
Cold, dry summer, very cold winter 
Cold, dry winter, hot summer 

Cold, dry winter, warm summer 

Cold, dry winter, cold summer 

Cold, dry and very cold winter 
Cold, no dry season, hot summer 

Cold, no dry season, warm summer 

Cold, no dry season, cold summer 

Code 
Cwc 
Dsa 
Dsb 

Dsc 
Dsd 
Dwa 
Dwb 

Dwc 

Dwd 
Dfa 
Dfb 

Dfc 

Cold, no dry season, very cold winter 

Dfd 

Page 37 of 52 

-23-23AfAmAwBShBSkBWhBWkCfaCfbCfcCsaCsbCwaCwbDfaDfbDfcDfdDsaDsbDscDsdDwaDwbDwcDwd123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100cfmnsbmnmxdmnpfmnslmnaxis1axis2-2.01.5-2.02.0TrAGiTrAGrTrBRmTrBRuTrFRmTrFRuTrGRmTrGRuTrSHmTrSHuTrSPmTrSPuTrWTL1234567891011121314151617181920212223242526272829303132333437384345485253545556575859606970717273757778808182838485868788899091939496979899100  SPECIES  SAMPLESis14 = 1is14 = 2is14 = 3is14 = 4is14 = 5is14 = 6is14 = 7is14 = 8is14 = 10is14 = 11is14 = 12is14 = 13is14 = 14 
 
 
 
 
 
 
 
 
 
 
 
 
Fit for the EFTs for 
the Köppen Zones 
(DCA)  

Figure 19: Percentage fit for the EFTs in the DCA with the Köppen zones. The values represent the % 
of  variance  in  all  Köppen  Zones  that  the  ordination  with  the  EFTs  explain.  Large  values  indicate 
good fit of the two datasets. 

In the DCA with the land use system classes the first four axes explained 64.3% of the 
variation in  the data, showing  good correspondence  with  the  ISODATA  clusters  (table 10). 
The first two axes explained 57.5% of the variance and as shown by the bi-plot (Figure 20) 
most of the ISODATA clusters could be associated to one of the land use classes (Figure 21). 
The CA axes were detrended by 2nd order polynomials as explained before.  

Table 10: Detrended Correspondence Analysis (DCA) of the ISODATA classes with the FAO land use 
system classes 

Eigenvalues 
Cumulative % of explained variance 

Axis1 
0.568 
39.4 

Axis2 
0.261 
57.5 

Axis3 
0.065 
62.0 

Axis4 
0.032 
64.3 

The above outlined characterization of  global ecosystems  is  based on an up-to-date, 
effective and repeatable indicator system obtained using Phenolo variables. The derivation of 
Ecosystem  Functional  Types  from  remote  sensing  data  offers  an  objective  and  repeatable 
method to characterize the functioning of ecosystems and provides a stratification that can be 
used to implement the LSN method.  

Page 38 of 52 

 
 
 
 
 
 
 
 
 
 
Figure 20: Detrended Correspondence Analysis biplot of the  ISODATA clusters classified in the 14 
dendogram  categories  (icons  and  colors)  and  the  FAO  land  use  system  classes  within  the  major 
climatic zones (crosses).  

Figure 21: Percentage fit for the EFTs for the DCA with the Koeppen zones. The values represent the 
% of variance in the FAO land use system classes that the ordination with the EFTs explain. Large 
values indicate good fit of the two datasets. (reprint from Ivits et al., 2013d) 

Page 39 of 52 

-1.52.0-1.52.0ArAGiArAGrArBRmArBRuArFRmArFRuArGRmArGRuArSHmArSHuArSPmArSPuArWTLCdAGiCdAGrCdBRmCdBRuCdFRmCdFRuCdGRmCdGRuCdSHmCdSHuCdSPmCdSPuCdWTLTmAGiTmAGrTmBRmTmBRuTmFRmTmFRuTmGRmTmGRuTmSHmTmSHuTmSPmTmSPuTmWTLTrAGiTrAGrTrBRmTrBRuTrFRmTrFRuTrGRmTrGRuTrSHmTrSHuTrSPmTrSPuTrWTL123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100axis1axis2-2.01.5-2.02.0TrAGiTrAGrTrBRmTrBRuTrFRmTrFRuTrGRmTrGRuTrSHmTrSHuTrSPmTrSPuTrWTL1234567891011121314151617181920212223242526272829303132333437384345485253545556575859606970717273757778808182838485868788899091939496979899100  SPECIES  SAMPLESis14 = 1is14 = 2is14 = 3is14 = 4is14 = 5is14 = 6is14 = 7is14 = 8is14 = 10is14 = 11is14 = 12is14 = 13is14 = 14 
 
 
 
 
2.4.2  Implementing the Local Net Scaling method  

The Local Net production Scaling (LNS) method estimates the potential or maximum 
production  in  homogeneous  land  capability  classes  using  remotely-sensed  observations  and 
then models the actual productivity of a pixel as deviation from this local potential (Prince et 
al.,  2009).  Variation  in  the  computed  potential  production  can  be  caused  by  differences  in 
land  use,  land  cover  and  physical  factors  whereas  the  variation  is  reduced  by  stratification 
into  homogeneous  regions.  The  actual  productivity  related  with  a  local  potential,  reflects  a 
current level of efficiency of productivity of an area and therefore this information is used as 
Productivity Status Map.  

The  homogeneous  units  were  computed  as  in  section  2.4.1  representing  global 
Ecosystem  Functional  Types  (EFTs)  derived  from  a  29  years  remote  sensing  time-series 
(1982-2010). The proportion between the potential and actual productivities was expressed in 
terms of the Cyclic Fraction  variable (Figure 3 and Figure 22)  rather than using the remote 
sensing  approximated  total  production  of  the  ecosystems  or  Standing  Biomass.  The  Cyclic 
Fraction  is  subject  to  high  inter-annual  variability  that  might  be  due  to  both  natural  and 
anthropogenic factors that are representative or conditioning the variability of this net annual 
productivity. The Cyclic Fraction expresses that part of total Standing Biomass productivity 
that  is  potentially  regularly  appropriated  by  human-environmental  use,  hence  it  is  a  good 
approximation  to  denote  current  productivity  levels.  In  order  to  spatially  accommodate  a 
representative number of pixels within the EFTs (Figure 23) the Cyclic Fraction was derived 
from the SPOT VEGETATION sensor on 1km spatial resolution. The yearly Cyclic Fraction 
values  were  computed  for  the  years  2006-2010  with  the  Phenolo  software  and  then 
temporally  averaged  into  one  value  over  this  5  year  period  (Figure  24).  This  5  year  period 
value  represents  the  average  annual  production  or  current  productivity  status  that  overlaps 
with  the  end  of  the  1982-2010  time  series  that  was  used  to  compile  the  Long  Term 
Productivity Change Map as explained above.  

Within each EFT the 90 percentile of the Cyclic Fraction values were computed and 
were  considered  representative  for  the  potential  or  maximum  level  of  annual  production  of 
that  given  homogeneous  spatial  unit.  Values  larger  than  the  90  percentile  were  considered 
outliers and were replaced by this 90 percentile value. Following the LNS approach, for all 
pixels within the EFT unit the proportion of the pixel value and the 90-percentile value was 
computed.  This  value  stands  for  the  local  scaled  productivity  of  that  pixel.  The  resulting 
values  are  between  0  and  1  where  values  close  to  1  represent  pixels  with  Cyclic  Fraction 
close to the local ecosystem’s highest annual production. In that sense for each pixel a level 
of  efficiency  as  compared  to  its  local  potential,  expressed  in  percentage  (listed  value  in 
Figure 25 times 100) performance is obtained. Figure 25 illustrates spatial results for Europe.  

This information as proxy for the current status of efficiency of productivity will be 
used jointly with the before explained long-term dynamics in view of compiling a base layer 
for further land degradation assessment.  

Page 40 of 52 

 
 
 
 
SB 

Figure  22:  Schematic  representation  of  Standing  Biomass  (are  under  the  curve)  and  the  Cyclic 
Fraction determined by the vegetation growing season.   

Ecosystem Functional Types (EFTs) 

Figure 23: Spatial units of the Ecosystem Functional Types (EFTs) of Europe 

Page 41 of 52 

 
 
 
 
 
Figure 24: Cyclic Fraction or annual productivity averaged for the years 2006-2010 calculated from 
SPOT Vegetation time-series data. 

.  

Page 42 of 52 

 
 
 
 
 
Figure 25: Local Net Scaling of the Cyclic Fraction within the Ecosystem Functional Types polygons. 

Page 43 of 52 

 
 
2.5  Combined assessment of land-productivity 

Analysis of long-term changes and current efficiency levels of vegetative Standing Biomass 
are combined into land-productivity dynamics according to the scheme presented in figure 1. 

Output from both the Long-Term Change Map, based on the Steadiness Index was combined 
with  start  levels  and  state  change  of  productivity  (Figure  15),  and  with  the  potential 
productivity map (Figure 25), based on the Local Net Scaling approach. This Local Scaling 
of the Cyclic Fraction was considered in two classes:   

1: less than 50 % of the EFTs highest Cyclic Fraction 
2: more than or equal to 50 % of the EFTs highest Cyclic Fraction 

Both the change and the status maps  are qualitative classification schemes. Combining these 
two together will therefore be based on a lookup table classification as shown in table 11. 

Figure 26 shows the final result of the combination of long-term change and status map  
based on the classification scheme of table 11. This final result shows the land-productivity 
dynamics over the EU.  Global maps are under production..  

Figure 26 maps five classes of land-productivity levels over the EU: 

Declining Land-productivity 
Early-signs of decline of land-productivity 
Stable, but stressed land-productivity 
Stable, not stresses land-productivity 
Increasing land-productivity 

For instance the class “declining land-productivity dynamics” is assigned to areas that  

(a) have showed evidence of prevailing downwards trends in standing biomass over 
the  twenty-nine  year  observation  period  1982-2010,  using  the  steadiness  index 
combined with  starting levels  and state change of productive levels (see Figure 15), 
and in combination show  
(b)  a  current  lower  than  potential  production  efficiency,  based  on  a  productivity 
compared  to  its  contextual  maximum  using  higher  resolution  data  for  a  baseline 
window of 5 years at the end of the time series (2006-2010).  

Hence,  land-productivity  dynamics  can  indicate  levels  of  sustained  land-quality  and  is 
therefore used as first step in the land degradation assessment. 

Page 44 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
Long Term Productivity Change Classes 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

St1 
lo 
- 

St1 
lo 
1 

St1 
lo 
2-> 

St1 
me 
- 

St1 
me 
1 

St1 
me 
2-> 

St1 
hi 
- 

St1 
hi 
1 

St1 
hi 
2-> 

St2 
lo 

St2 
me 

St2 
hi 

St3 
lo 

S32 
me 

St3 
hi 

St4 
lo 
- 

St4 
lo 
1 

St4 
lo 
2-> 

St4 
me 
- 

St4 
me 
1 

St4 
me 
2-> 

St4 
hi 
- 

LS  < 
50% 

LS  >= 
50% 

1 
d 

2 
ew 

1 
d 

2 
ew 

1 
d 

1 
d 

1 
d 

2 
ew 

1 
d 

2 
ew 

1 
d 

1 
d 

2 
ew 

3 
nf 

1 
d 

2 
ew 

1 
d 

2 
1? 

3 
nf 

3 
nf 

3 
nf 

3 
nf 

3 
nf 

3 
nf 

4 
pf 

4 
pf 

4 
pf 

4 
pf 

4 
pf 

4 
pf 

4 
5 

5 
i 

4 
5 

5 
i 

5 
i 

6 
si 

4 
5? 

5 
i 

5 
i 

5 
i 

5 
i 

6 
si 

5 
i 

6 
si 

Column headers (see figure 15):  
St1-4: the four steadiness classes 
lo: low; me: medium; hi: high Standing Biomass productivity at the beginning of the time-series  
-: no class change; 1: changed 1 class; 2: changed 2 or more classes 
Table values: 
1-5: the five classes of land productivity levels 
d:  declining  land  productivity;  ew:  early-signs  of  decline  of  land  productivity;  nf:  negative  fluctuation  (stable,  but  stressed  land-productivity);  pf:  positive 
fluctuation (stable, not stressed land-productivity); i:increase (increasing land productivity); si: strong increase (increasing land producitivity) 
Row headers: 
LS = Local Scaling; values are as in Figure 25.  

Table 11: Look up table for the combination of the combined steadiness and the local scaling of the Cyclic Fraction 

Page 45 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figure 26: An assessment of land productivity dynamics in Europe combining productivity 
dynamics of Standing Biomass and the local scaling of the Cycling Fraction. 

Page 46 of 52 

 
 
 
 
 
 
3  Conclusion and Outlook 

Land  degradation  and  desertification  assessment  are  very  complex  and  there  is  no  agreed 
scientific methodology or measurement protocol  today that  is  applicable  at the  global level. 
Satellite observations offer potential to monitor that part of the earth’s cover that is or can be 
used as ecosystem services responding to human demands.  

This  report  documents  how  land-productivity  dynamics  can  be  calculated  from  vegetation 
indices  derived  from  long-term  low-resolution  satellite  time  series  such  as  the  GIMMS3g 
dataset  combined  with  productivity  efficiency  measurements  derived  from  short,  recent, 
medium resolution data such as those from the SPOT VEGETATION sensor.  

Land  degradation  oriented  studies  before  have  been  based  on  NDVI  measurements  mainly. 
Here  this  limited  set  of  variables  is  extended  by  de-convolution  of  the  original  times  series 
data  resulting  in  a  number  of  phenological  and  productivity  relevant  variables.  Keeping 
statistical  solidity,  qualitative  approaches  have  been  selected  to  classify,  interpret  and 
integrate  these  several  variables.  Long-term  change  dynamics,  expressed  by  the  Steadiness 
index, have been combined with productivity levels at the origin of the time series, and with 
information on state change. This long-term dynamics are then combined with information on 
the current status of efficiency of productivity. This calculation has been done considering the 
local contextual potential or maximum levels of productivity and gives an impression of how 
good or bad a certain land-area is functioning.  

Considering  the  data  resolution  and  combined  processing  as  explained,  it  is  felt  that  land-
productivity dynamics as a holistic approach reflects climatic constraints, the overall quality 
of the land, the efficiency of using the soil and other resources and it indirectly indicates the 
level  at  which  these  resources  are  appropriated  for  human  use,  i.e.  is  the  land  used  for 
intensive/extensive  agriculture,  for  grazing,  forestry  or  urban?  Hence,  land-productivity 
dynamics can indicate levels of sustained land-quality and is therefore used as first step in the 
land  degradation  assessment.  Large-scale  land-productivity  dynamics  is  good  base  level 
information  on  which  to  integrate  further  contextual  knowledge  for  finally  assessing  land 
degradation  and  desertification  conditions.  This  integration  analysis  is  needed  to  obtain  a 
holistic  interpretation  of  possible  on-going  land  degradation  that  comprehends  and  explains 
the biophysical dynamics in relation to anthropogenic aspects. 

Further  work  will  focus  on  methods  for  integrating  ancillary  data.  A  first  example  of  such 
exercise is given in the Annex where from the base land-productivity dynamics layer further 
analysis are implemented to understand, explain and interpret results in view of  contributing 
to land degradation assessment frameworks.  

Using and interpreting these data, a comprehensive report on “land-Productivity Dynamics in 
Europe - Towards a Valuation of Land Degradation in the EU” has been compiled. 

Page 47 of 52 

 
 
 
 
 
 
 
 
Page 48 of 52 

 
 
References 

Bradley, R. Jacob, J. Hermance, J.F. Mustard, 2007. A curve-fitting technique to derive inter-
annual  phonologies  from time series  of noisy  NDVI satellite data.  Remote Sens.  Env., 106, 
pp. 137-145. 

FAO. Land Degradation Assessment in Drylands: Mapping Land Use Systems at Global and 
Regional  Scales  for  Land  Degradation  Assessment  Analysis  v1.1.,    Food  and  Agriculture 
Organization of the United Nations: Rome, Italy, 2011. 

Fensholt, R., Proud, S. R. Evaluation of Earth Observation based global long term vegetation 
trends  — Comparing GIMMS and MODIS  global  NDVI time series.  Remote Sens Environ. 
2012, 119, 131–147. 

Guo, W.Q., Yang, T.B., Dai, J.G., Shi, L., Lu, Z.Y., 2008. Vegetation cover changes and their 
relationship to climate variation in the source region of the Yellow River, China, 1990-2000. 
International Journal of Remote Sensing, Vol. 29, No. 7, 2085-2103. 

Hill,  M.O.,  and  Gauch,  H.G..  Detrended  correspondence  analysis,  an  improved  ordination 
technique. Vegetation 1980, 42, 47-58. 

Ivits, E., Cherlet, M., Sommer, S., Mehl, W. 2012a. Ecosystem Functional Units characterized 
by  satellite  observed  phenology  and  productivity  gradients:  a  case  study  for  Europe. 
Ecological Indicators 27, 17–28. 

Ivits, E., Cherlet, M., Toth, G., Sommer, S., Mehl, W., Vogt, J., Micale, F., 2012b. Combining 
satellite  derived  phenology  with  climate  data  for  climate  change  impact  assessment.  Global 
and Planetary Change 88–89, 85–97. 

Ivits, E., Cherlet, M., Sommer, S., Mehl, W. 2012c. Addressing the complexity in non-linear 
evolution  of  vegetation  phenological  change  with  time-series  of  remote  sensing  images. 
Ecological Indicators 26, 49–60. 

Ivits,  E.,  Cherlet,  M.,  Horion,  S.,  Fensholt,  R.,  2013d.  Global  Biogeographical  Pattern  of 
Ecosystem Functional Types for Land-Degradation Studies. Remote Sensing, XXXX 

Jönsson,  P.  and  Eklundh,  L.,  2004.  TIMESAT-A  program  for  analyzing  time-series  of 
satellite sensor data. Computers & Geosciences, 30, pp. 833-845 

Linderholm, H.W., 2006. Growing season changes in the last century. Agricultural and Forest 
Meteorology, 137, 1-14. 

Myneni, R.B., Keeling, C.D., Tucker, C.J., Asrar, G. and Nemani, R.R. 1997: Increased plant 
growth in the northern high latitudes from 1981 to 1991. Nature, 386: 698-702. 

Parmesan, C., 2006. Ecological and evolutionary responses to recent climate change. Annual 
Review of Ecology Evolution and Systematics, 37,637-669. 

Paruelo, J.M., Jobbagy, E.G., Sala, O.E. Current distribution of ecosystem functional types in 
temperate South America, Ecosystems 2001, 4, 683-698. 

Page 49 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
Peel,  M.C.,  Finlayson,  B.L.,  McMahon,  T.A.  Updated  world  map  of  the  Köppen-Geiger 
climate classification. Hydrol. Earth Syst. Sci. 2007, 11, 1633-1644. 

Prince, S.D. Mapping desertification in Southern Africa. In: Gutman, G., Janetso, A., Justice, 
C.O., Moran, E.F., Mustard, J.F., Rindfuss, R., Skole, D., Turner, B.L. (Eds.). Land Change 
Science:  Observing,  Monitoring  and  Understanding  Trajectories  of  Change  on  the  Earth's 
Surface. Springer, Berlin, 2004, pp. 163-184. 

Prince,  S.D.,  Becker-Reshef,  I.,  Rishmawi,  K.,  2009.  Detection  and  mapping  of  long-term 
land  degradation  using  local  net  production  scaling:  Application  to  Zimbabwe.  Remote 
Sensing of Environment 113 (2009) 1046–1057. 

Reed, B. C., Brown, J.F., VanderZee, D., Loveland, T. R., Merchant, J. W.,
Ohlen,  1994,  Measuring  phenological  variability  from  satellite
Vegetation Science, 5: 703-714. 

and D. O 
imagery,  Journal  of 

Sommer, S.,  C. Zucca, A. Grainger, M. Cherlet, R. Zougmore, Y. Sokona, J. Hill, R. Della 
Peruta,  J.  Roehrig,  G.  Wang,  2011:  Application  of  indicator  systems  for  monitoring  and 
assessment of desertification from national to global scales. Land Degrad. Develop. 22: 184–
197 (2011).TerBraak, C.J.F. and Prentice, I.C. A theory of gradient analysis. Adv. Ecol. Res. 
1988, 18, 93-138. 

Stow,  D.;  Daeschner,  S.;  Boynton,  W.;  Hope,  A.  Arctic  tundra  functional  types  by 
classification of single-date and AVHRR bi-weekly NDVI composite datasets. Int. J. Remote 
Sens. 2000, 21, 1773–1779. 

TerBraak,  C.J.F.  and  Prentice,  I.C.  A  theory  of  gradient  analysis.  Adv.  Ecol.  Res.  1988,  18, 
93-138. 

Wessels,  K.J.,  prince,  S.D.,  Malherbe,  J.,  Small,  J.,  Frost,  P.E.,  VanZyl,  D.  Can  human-
induced  land  degradation  be  distinguished  from  the  effects  of  rainfall  variability?  A  case 
study in South Africa. J. Arid Environ. 2007, 68, 271-297. 

White, M.A., DeBeurs, K., Didan, K., Innouyes, D., Richardson, A. D., Jensen, O., O’Keefe, 
J., Zhang, G., Nemani, R. R., VanLeeuwen, W. J. D., Brown, J. F., DeWit, A., Schaepman, 
M., Lin, X., Dettinger, M., Bailey, A.S., Kimball, J., Schwartz, M.D., Baldocchi, D. D., Lee, 
J.  T.,  Lauenroth,  W.  K.,  2009.  Intercomparison,  interpretation  and  assessment  of  spring 
phenology in North America estimated from remote sensing from 1982-2006. Global Change 
Biology, 15, 2335-2359. 

X.  Zhang,  M.A.  Friedl,  C.B.  Schaaf,  A.H.  Strahler,  J.C.F.  Hodges,  F.  Gao  et  al.,  2003. 
Monitoring vegetation phenology using MODIS. Remote Sens. Env., 84, pp. 471-475. 

Page 50 of 52 

 
 
 
 
 
 
 
 
 
 
 
 
 
Europe Direct is a service to help you find answers to your questions about the European Union 
Free phone number (*): 00 800 6 7 8 9 10 11 
(*) Certain mobile telephone operators do not allow access to 00 800 numbers or these calls may be billed. 

A great deal of additional information on the European Union is available on the Internet. 
It can be accessed through the Europa server http://europa.eu 

How to obtain EU publications 

Our publications are available from EU Bookshop (http://bookshop.europa.eu), 
where you can place an order with the sales agent of your choice. 

The Publications Office has a worldwide network of sales agents. 
You can obtain their contact details by sending a fax to (352) 29 29-42758. 

Page 51 of 52 

 
 
 
 
 
 
 
 
-

L
B
N
A

-
2
6
0
5
2
-

E
N
N

-

JRC Mission 

As the Commission’s  
in-house science service,  
the Joint Research Centre’s  
mission is to provide EU  
policies with independent,  
evidence-based scientific  
and technical support  
throughout the whole  
policy cycle. 

Working in close  
cooperation with policy  
Directorates-General,  
the JRC addresses key  
societal challenges while  
stimulating innovation  
through developing  
new methods, tools  
and standards, and sharing  
its know-how with  
the Member States,  
the scientific community  
and international partners. 

Serving society  
Stimulating innovation  
Supporting legislation 

Page 52 of 52