Today we will have a Webinar that presents 4 topics.
As a fifth topic, your Questions as typed in the Q&A-box, will be addressed.
First, the basics of mapping Agro-EcoSystems will be presented,
after which I will show the strengths of using the dimension ‘time’ for mapping and monitoring Agro-EcoSystems.
Then, a series of maps derived from NDVI data-cubes will be presented,
followed by more details on data sources, pre-processing, visualization, and interpretation.
Section-1. Hyper-Temporal RS gives you a rich data-cube that can be interpreted
to provide information on what is growing where,
when it is growing, how it is performing, and when/how it is changing.
This information is essential: ...
The Index
NDVI is the most widely used and possibly most successful ‘spectral difference’ index.
Which is the ratio of: the difference of the near infrared and red reflection, over their sum.
We see in the left figure that the green line, which represents vegetation
[factually, it represents the amount of active chlorophyl present],
absorbs most light in the visual red band and reflects most light in the Near Infra Red [NIR] band.
By comparing both bands, we can thus measure the amount of active chlorophyl.
An NDVI data-cube
A data-cube can be presented as a movie.
The long-term median statistics (or: climatology) of actual (cleaned) NDVI-data are presented in the left figure.
The movie consists of 36 sequential images,
each representing a specific dekad, which is a 10-day period.
The figure displays a high and dense amount of vegetation in orange/red,
and bare-soil (no vegetation) as dark-blue/black.
In between light to dark-green colours are used.
The movie clearly shows spatial-temporal changes in vegetation phenology.
These are caused by climate differences, soil/terrain/landform differences,
land use differences, etc.
All differences are reflected through land cover (read: NDVI) differences.
Because of these differences, a specific land use, for instance:
cropping of common wheat, can be mapped.
In the right figure, an example of such a map is presented.
The “common-wheat crop-intensity map” is produced on the basis of this movie,
on available agricultural statistics (by administrative areas),
on field-observations (the Copernicus-Lucas data),
and the EU-Corine map (based on a 30m TM-imagery interpretation).
Thus, satellite imagery and tabular data can be fully integrated and transformed into maps
containing information of importance.
Let us recapture:
Land Cover (NDVI) relates, amongst others, to:
Since many of the above have a strong temporal behaviour,
also the NDVI varies over time.
Often, clear crop growing periods (seasonality) can be recognized
through plotting NDVI over time; see the top-right figure,
which displays a specific NDVI-Profile.
Such NDVI-profiles may differ over years because crop-performance will differ
due to the severity of adverse growing conditions (=perils).
The bottom-right figure displays such annual (seasonal) differences.
Extreme heat, combined with drought, caused severe crop-failure during the 2010 season;
2003 however delivered a bumper crop.
If perils hit a relatively large agricultural area, it often leads to severe price hikes followed by unrest.
Section-2 will zoom-in on differences in behaviour of individual pixels,
but also of specific areas, over time.
Comparing NDVI-profiles between pixels, between areas, and between seasons (years)
provides useful information on the actual land cover present,
of the variability of the vegetation performance,
and of gradual to abrupt changes in land cover.
Visual inspection of NDVI time-series by pixel
The shown visual inspection of 3 individual pixels clearly shows that land cover of the 3 pixels truly differs.
The shown NDVI-Profiles (DN-values) cover the period 2000-2013
(13 years of MODIS-Terra imagery).
- The red-line displays a continuously high NDVI-amount (DN-values) that indicates an evergreen situation.
Factually, the pixel is located in the Mau Forest National reserve of Kenya,
which consists predominantly of a dense forest complex which includes dense bamboo patches.
- The dark-green line displays clear seasonality with 2 peaks annually indicating 2 crop-growing seasons.
The respective landscape is dominated by arable fields with dispersed some trees
- The light-green line displays a 1-season per year NDVI-profile,
with relatively low NDVI-values indicating sparse vegetation.
On the respective Google-Earth image, a bare-landscape is shown;
likely, seasonal vegetation (short grasses) dominate the specific area.
Each pixel will have its own ‘history’
…and, each pixel will have its own ‘variability’
Each area will have its own ‘variability’
Plotting on top of each-other all 13 annual NDVI-Profiles (DN-values),
representing whole areas in and around the Masai Mara Park in Kenya,
displays how reliable the seasonality of the respective areas are.
Inspect: - If there are shifts in timing of seasons causing early or late season starts and/or season ends, and
Inspect: - If complete absence of seasons do occur.
Where area 48, in general, is dryer and with less vegetation than area 35,
both areas only confirm that from April to June a clear ‘wet’ season occurs.
The remainder of the year can remain completely dry or can experience another early or late wet season (with every option in-between).
The shown unreliability in seasonality makes for this area arable agriculture rather difficult;
good years will be fully mixed with years experiencing severe crop failures.
Semi-extensive animal husbandry with local trans-migration options is thus the best option to build ones livelihood.
The Masai thus got it right !
Monitoring an area through its NDVI data
… The drop in cultivated area was caused by a major overhaul of the irrigation infrastructures upstream.
Thus, the shown ‘drop’ was only temporary.
Monitoring an area through its NDVI data
Another example.
Section-3. The third section shows several mapping-examples, and reviews a series of aspects that relate to creating area-strata from NDVI data-cubes.
The mapping examples cover: ...
Pixels are Mixels
This example relates to wheat & barley systems in Andalusia, Spain.
When an NDVI data-cube is converted to a strata map,
and if only the key-strata with wheat and barley areas are extracted,
then we get the shown 7 NDVI-profiles.
It shows profiles of 5 wheat strata and of 2 barley strata.
Where barley profiles mainly differ due to the cropping intensity (22 versus 13% in average),
the wheat strata differ both in wheat cropping intensities (29 to 47%),
and in the presence of sunflower cropping (22 to 40%).
As an inset, a map of Andalusia, showing the overall wheat cropping intensity (area-fractions) is presented;
it varies, if present, from 2 to 47%.
A key message of this slide is that preparing a ‘library’ of NDVI-profiles makes no sense.
Crop mixes vary, responses over time vary, etc. etc.
A profile will be very-very area-and-season specific and very specific in what mixes it represents.
We need thus skills to stratify a space-time cube and to relate each created strata with secondary data
(or as a base to carry out fieldwork)
in order to make a legend as required.
Mapping at which resolution ?
A general misconception is that use of higher spatial resolution data, produces better mapping opportunities.
However, please note that the sought-after strata must differentiate land cover & land use patterns
in the space-time domain that relates best to differentiating different climate, soil, terrain, and landforms.
Differentiating dominant mixtures of Agro-EcoSystems happens to relate better to landscapes when using coarse resolution pixels.
‘Mixels’ at, for instance 100m resolution, will be less ‘generalised’ and will produce an increasing mix of land cover type combinations.
Note also that we mostly analyse NDVI-based data-cubes at country and regional scales.
Using a too fine spatial resolution of such huge areas will also defeat the effort to create generalizations at the required scale.
As example, for the Mekong area of Vietnam, two stratified maps are presented,
the left one at 1km resolution (using Spot-VGT)
and the right one at 250m resolution (using MODIS-Terra.
The 1km map does not suffer from various inclusions (specific mixels) within mapped landscapes
and provides a better generalization (as required).
The right map shows a serious increase in specific mixels,
each with a different unique combination of land cover types.
As such, linear features like canals and roads, create many of these mixels.
They clearly interfere with the generalization effort envisaged.
Mapping the gradient between a desert and continental humid areas
Stratifying, for instance an NDVI data-cube of Mongolia, creates a series of strata that follow a beautiful gradient:
from desert and steppe to continental humid (severe cold winters and warm/hot summers).
Capturing such gradients is simply impossible without the use of an NDVI-based spatial-temporal data-cube.
A down-hill gradient
Another nice example to map a gradient occurs in Java, Indonesia.
The white arrow points to an area East of Jakarta.
In the whole area cropping occurs during two seasons annually.
The timing of each season gets delayed when one moves to the North, or better: when one moves downhill towards the Java Sea.
This specific gradient is divided into 4 strata, but it could have been divided into more strata,
each with a ‘narrower’ temporal window regarding the start and end dates of each season.
Exactly, because of the existence of many criss-crossing gradients all over the globe,
it becomes impossible to specify “how many strata to differentiate”.
It is as difficult as to state that there just x-shades of grey (between white and black).
Just consider that you aim to create a map (for monitoring) which must contain a certain level of ‘generalization’.
Only the expert (you) can decide on the required level (pending on your objectives).
Base-Line Mapping
Base-line mapping, through an NDVI-based space-time data-cube,
at country or regional level,
is meant to differentiate areas that are clearly different regarding land use and land cover.
The created strata or map-units must be relatively homogeneous within regarding the land cover mixes present,
and rather heterogeneous between the strata.
A base-line map can be used, amongst others, for:
- To adopt the strata as a first-level survey sample scheme (i.e. stratified sampling),
e.g. for Area Frame Sampling (AFS) to collect data that relate to land cover (biodiversity) and/or to land use (agricultural statistics),
- For developing a subsequent monitoring scheme to detect changes in land cover and land use,
- For subsequent mapping individual strata at higher spatial scales through use of high resolution earth observation imagery,
- To use the created NDVI-profiles by strata for subsequent monitoring the performance of those strata as required e.g. for early warning or for insurance schemes
Note that mapped ”forests” of Korea (top-right figure with the green lines)
differentiates a clear gradient capturing differences in forest-densities (as average for a 1km pixel).
Note also that the ‘dent’ in the green NDVI-profiles reflects the fact that during the ‘long’ growing season of trees,
haze and clouds influenced (dented) the NDVI-readings.
The ‘dent’ is thus an artifact specific for the forests and the prevailing weather conditions.
A second example of a base-line map concerns Pakistan.
Using 11 years of MODIS NDVI-imagery at 250m resolution,
80 different clusters were created.
Subsequently, the NDVI-profiles of the 80 clusters were grouped into 13 major-groups.
The map shows that the colour scheme of the Agro-Ecosystems map of Pakistan follows the 13 colours used for plotting the NDVI-profiles (legend).
This map clearly reflects (i) an altitude gradient (up North),
(ii) Irrigated areas versus Barani (rainfed) areas,
(iii) the performance of the two cropping seasons per year.
From ‘purple’ to ‘red’, to ‘blue’, the performance of the first season was best,
then both seasons performed equally,
and more South (downstream) the 2nd season performed best.
The 3rd Base-Line map concerns a specific land use:
rice-cropping in India and Pakistan.
The rice-systems are differentiated based on
(i) the water management system, and
(ii) differences in crop calendars practiced.
The latter differences were nicely captured through the respective area-specific NDVI-profiles.
A small selection of these profiles is shown;
note specifically the ‘sudden dip’ in NDVI at the start of a growing season
as caused by flooding of the concerned paddies.
The last base-line mapping example concerns an effort by IWMI-staff, Colombo,
to map the global cropping areas of five major crops.
They are: wheat, rice, maize, barley and soybeans.
Also in this case, the authors made ample use of differentiating NDVI-profiles.
Using the link provided at the bottom of the map will play a video by ESA titled “Changing Lands”.
This video further illustrates the strength to map actual land cover and to monitor actual land cover
through the use of remotely sensed imagery ‘over time’.
Topic-4. Our 4th topic will present a series of examples to illustrate:
acquisition, pre-processing, visualization, analysis and interpretation,
of Hyper-Temporal NDVI-based data-cubes.
But first, let us list some typical characteristics of Hyper-Temporal (HT) NDVI-catalogues.
They are: ...
Maximum Value Composites (=MVCs)
MVC means: The ‘best’ pixel-specific, daily, NDVI-value of a given period is used
to represent that specific period.
This follows the logic that haze and clouds ‘allways’ reduce the actual NDVI-value
which leads to a lower ‘recorded’ NDVI-value.
As it happens, the highest recorded NDVI-value in a given period
thus suffered least from haze and cloud contamination (of course assuming that during the given period,
the NDVI remained relatively static).
Since during a 10-day period, the ‘greenness’ of vegetation is not supposed to fluctuate heavily,
use of a 10-day MVC-period is very common.
16-day and monthly composites already (partially) defeat the purpose of creating a Hyper-Temporal NDVI-based data-cube.
The very best 10-daily MVC NDVI-series, as made available by the Copernicus Global Land Service (CGLS), are:
For the period 1 Jan.’1999 – 30 June 2020: The 1km SPOT-VGT & Proba-V catalogue (version 3; final version!), and
For the period 1 July 2020 – present: The 300m Sentinel-3 OLCI-based catalogue (version 2).
Both series are fully cross-calibrated (BRDF-adjusted), and carry a guarantee that they will be continued till at least 2040.
Why a BRDF-adjusted NDVI ?
The overpass times of the various platforms used, differed and varied.
See for instance the ‘drifting’ of Proba-V in the figure on: “The evolution of local overpass time at nadir”.
In addition, the viewing angle of the overpassing satellite to the specific ‘pixel’ varies between image rows & columns,
and pending calendar dates, to make it worse, also the sun follows an annual North-South cycle.
This all impacts on absorption-reflection aspects of VIS and NIR light.
It also impacts on the shade pattern thrown by objects ‘catching’ the sunlight.
To make it worse, also the slope and slope-position of the surveyed ‘pixel’ matters.
All the above is adjusted by the data-provider using a Bidirectional_Reflectance_Distribution_Function (BRDF).
The bottom-right figure shows the results of the BRDF exercise, while the top-right figure shows the considerable mismatch between catalogues when a BRDF-adjustments is omitted.
… more pre-processing by ‘you'
First note that in the CGLS-imagery, all DN-values above 250 are Quality Flags,
and when used for a pixel, that pixel has thus a ‘missing’ NDVI-value.
Then after the required BRDF-adjustment, and after creating the MVC,
the resulting NDVI-data can still be influenced by ‘rather persistent’ haze and clouds.
The brown-line in the shown figure, displays for a specific pixel,
the NDVI data-series covering 2006.
Many negative ‘outlyers’ and downward anomalies do occur.
Theoretically, the NDVI pattern, covering the specific pixel and year, would follow the dotted green line.
That line can partially be estimated through a function called ‘Timesat’ (by Ecklund et.al).
Timesat is iterative and uses the Savitsky-Golay smoothing algorithm to interpolate between NDVI-readings.
If the estimated ‘smoothing’ value is higher than the actual value at time-x, then Timesat uses the higher value.
This process is repeated thrice with a growing temporal window.
It results in an estimation of the “Upper-Envelope”.
At the peak of the growing season, the red-line did not fully manage to approach the dotted green line;
persistent haze and clouds made that impossible, and the Upper-Envelope is thus never completely perfect.
The TIMESAT filter will deliver a ‘clean’ data-cube
Besides filling values for missing records and adjusting values for negative outliers and anomalies,
Timesat will also remove smaller ‘peaks’; see the figure.
This removal is achieved through a final iteration,
in which the actual estimate of the last Savitsky-Golay smoothing is used and not the possible ‘higher’ actual.
Real “Upward” NDVI-anomalies occur only seldomly.
We call it ‘Glint’.
It happens when the distance between the pixel and the sensor is very high and has a very low angle
(the ‘signal’ has to pass lots of atmosphere);
in addition, it only occurs when the pixel consists of snow or ice.
This type of anomaly thus occurs rather rarely, and should not concern us the least.
Once a BRDF-adjusted data-cube of MVC-values has passed the Timesat filter, we can start using that data-cube for our ‘science’.
Once cleaned, do visualize the NDVI-imagery as a movie
Viewing an original cleaned data-cube is of importance, because ‘data-implosion’ through clustering or classification is not yet done.
Viewing the data-cube as a movie will display and show lots of specifics in which one might be interested.
Do take time to detect and understand features that later must also show in the analysis results.
The movie can be ‘sharpened’ by overlaying a semi-transparent hill-shade image (prepared from a 30m resolution SRTM).
Also add roads, rivers, etc. as required to identify and recognize features of the movie.
Shown is a ‘movie’ of Mozambique prepared through 36 images representing 36 sequential dekads (these are 10-day periods).
The NDVI-value of each pixel of a specific dekad-image represents the median value of all annual repeats of the cleaned NDVI-values for that dekad and pixel.
The shown ‘movie’ thus represents the ‘climatology’ of Mozambique:
it is a cloudless representative of the average green vegetation present over time covering the past 20 years (2000-2019).
It is thus the ‘best’ base-line country-level map to study present differences and monitor future changes in land cover and land use.
After one has studied the movie, it is time for a computer-guided ‘data-implosion’.
The ‘movie’ must become a 2D-map containing map-units, and its legend must become a collection of NDVI-profiles.
Each NDVI-profile represents one specific map-unit of the prepared map.
Note that, through the clustering or classification process of a time-series of data by pixel,
the pixel has received a map-unit code and the pixel-data (its time-series of NDVI-values)
has shifted to the legend.
Actually: each NDVI-profile in the legend represents the averaged values of all pixels classified to that specific cluster or class.
The produced map is an intermediate-map, ready for use during fieldwork or for integration with secondary data.
Note that this map can, for 100%, be reproduced at will.
The only influence one had is by defining the number of clusters that the 2D map must contain.
That represents the only ‘unknown’ when clustering (classifying) the movie.
Supervised classification, based a list of a-priori specified land-cover types,
is rather impossible to carry out, because the ‘coarse’ pixels will often contain complexes of land-cover types,
and one may easily overlook unique land-cover types that an unsupervised clustering algorithm like ISODATA will properly identify.
Note in the legend the “gradient” of NDVI-profiles.
Their mutual distances indicates that preparing 50 clusters ‘did the job’.
Another ‘movie'
Shown now is the climatology of Luzon, The Philippines,
representing the period 1998-2011.
Using the Iterative Self-Organizing Data Analysis Technique (ISODATA),
an unsupervised clustering of the Hyper-Temporal data-cube is performed.
Using patterns identification and classification, the ‘movie’ is ‘imploded’ to a 2D-map with 79 map-units;
each map-unit is described through an NDVI-profile.
Shown now are some NDVI-profiles that clearly differ and that can each be labelled using some expert knowledge and field-experience in Luzon.
The red line, for instance, shows 2 distinct growing periods within one year;
they suggests arable-cropping.
Indeed, in Luzon, the ‘red’ areas happen to be cropped to irrigated rice twice annually.
The light- and darker-green profiles show a distinct “dry-period”
and a rather long (extended) growing season, during which NDVI values are rather high.
Such findings relate likely to evergreen cover-types with a high density of trees and bushes.
Likely, the orange-line represents something in-between the red and green lines.
It depicts just one growing season. Likely it represents a mixel of rainfed rice and trees and/or bushes.
Verification through high resolution imagery and fieldwork is clearly required.
Gradually, through use of field data, expert knowledge, and secondary data,
the legend of the intermediate map will become a legend that contains factual land cover and land use terminology.
The intermediate NDVI-profiles will thus disappear.
Detecting past changes
In case, annual repeats are NOT generalized using their median values, and thus the movie does not represent ‘climatology’,
ISODATA clustering will produce a legend with NDVI-profiles covering the full time-period of the original data-cube.
Now, the variability between years is thus NOT removed from the dataset
- This also represents a choice made by the researcher
- Detecting ‘past’ changes or defining ‘past extremes’ will now be possible through studying the produced legend.
As example, shown are NDVI-profiles covering a 10-year period (1998-2008) and representing 3 specific irrigated rice areas in the Mekong area of Vietnam.
- In the blue area the extent of annual flooding in the area increased between 1998 and 2000.
Afterwards, a very steady sequence ensued of flooding and cropping twice a year.
- In the green area, the sequence of cropping and flooding was rather haphazard till early 2002.
Afterwards, annually, 3 crops were cultivated, with no time remaining for any serious flooding. Till 2007!
Then, famers claimed to have opened the dike that was completed in 2002 and that surrounds their villages,
to ‘harvest’ silt and thus to fertilize their paddies in a very natural way.
In was by intent and not an unforeseen disaster that created a loss of the specific crop.
- The brown line displays a very gradual increase in performance of the rice crop until 2003.
Otherwise, no shocks or abrupt changes did occur.
The last example will be presented through a sequence of slides.
It concerns mapping rice systems and their cropping calendars, in Odisha, India.
Step-1 shows the climatology movie, which was converted through ISODATA into a map with 100 clusters.
Step-2. Knowing that 100 clusters is rather ‘many’, the 100 NDVI-profiles were grouped into 21 distinct types of profiles.
This was partly achieved through k-means grouping using the profile-data,
and partly through ‘expert-interpretation’ considering both the ‘relevant’ parts of the NDVI-Profiles and the locations of related ‘profiles’,
that is: ‘where’ does a profile occur (and its extent).
This in-between slide shows that ‘expert-guided’ step by showing the spatial relation
of 4 created NDVI-groups, each having 2 distinct growing seasons:
- Groups [E-E1-I] occur closely together and they form together a gradient (see the profiles and the map)
- Group J (the red area) stands spatially separate and seems unrelated to the others.
Step-3 shows the NDVI-profiles of the 21 groups.
These profiles are used to identify the critical periods within one year that ISODATA used to differentiate the clusters (groups).
Inspection identified 4 distinct periods.
Knowing the 4 periods, a 2nd phase of ‘mapping’ started.
Steps 4 to 6
Step 4. For the identified 4-periods, NDVI-averages from all available cloud-free medium-resolution imagery,
covering recent years, were extracted.
This delivered 4 images that contained mean NDVI-values.
These 4 images were then used by ISODATA to once more generate 21 classes,
but now at a 30m spatial resolution.
The super-imposed figure illustrates some results of Step-4.
The 10 profiles shown all have just 4 datapoints each, matching the 4 earlier identified periods.
- Classes 11 & 15 relate to cropping during 2 seasons, where 11 occurs in the lowlands and 15 in the uplands.
- Classes 4, 6, 9 and 10 relate to cropping in the upland areas.
- Classes 5, 7, 8, and 12 relate to cropping in the lowlands.
Step 5: The extent of each of the 21 classes was inspected and compared with published seasonal area statistics on rice by administrative area (Ref. EARAS).
This created an improved legend containing area fractions by NDVI-class.
The shown table contains these area fractions as area-percentages, and does so by season, by water management system (irrigated versus rainfed), and by terrain unit (uplands versus lowlands).
Besides these fractions, also the fractions of the total net area sown annually to any crop is reported.
Please do note that in Odisha, used season-identifiers relate to the harvesting periods of the respective crops.
Step 6: The legend was further improved by matching created legend-items with published crop calendar information.
Note once more that the season-identifiers relate to the harvesting periods of the respective crops:
- Summer rice is harvested in June and has a Crop Growing Period Length of 120 to 135 days,
- Autumn rice is harvested in September, with a relatively short growing period of 80 to 100 days, and
- Winter rice is harvested from October to December, with a growing period that can extend from 100 to 150 days.
Finally do note, that I as author of this map, never set foot in Odisha.
Agronomic expert knowledge, great secondary data sources, and smart combination of NDVI-data in the spatial and temporal domains, did the job !
Section-5. This brings us to the end of the presentation.
Now it is time to address Questions posed by you in the Q&A chat-box.
As a repeat, the topics covered are repeated on this slide.