See discussions, stats, and author profiles for this publication at: [Link]

net/publication/238494196

Remote Sensing: Fundamentals, Types and Monitoring Applications of

Environmental Consequences of War

Article in Handbook of Environmental Chemistry, Volume 3: Anthropogenic Compounds · January 2009

DOI: 10.1007/978-3-540-87963-3_3

CITATIONS READS

4 8,902

2 authors:

Dhari al-ajmi Saif Uddin

Sultan Qaboos University Kuwait Institute for Scientific Research
38 PUBLICATIONS 280 CITATIONS 131 PUBLICATIONS 1,052 CITATIONS

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Climate Change Impacts in Gulf View project

Microplastics in Environment View project

All content following this page was uploaded by Saif Uddin on 21 October 2020.

The user has requested enhancement of the downloaded file.

Hdb Env Chem (2009): 41–123
DOI:10.1007/698_2008_6 41

Remote Sensing: Fundamentals, Types

and Monitoring Applications of Environmental
Consequences of War

Dhari Al Ajmi and Saif ud din

Abstract The chapter deals with the fundamentals of remote sensing, basic
principle of electromagnetic radiation and its interaction with the earth, atmosphere
and surface materials. The types of sensors, digital data formats, basic image
processing techniques, including image enhancement techniques and classification
methods are explained in details. Besides the image processing techniques, applica-
tions of remote sensing in Kuwait are presented in a section on applications, that
include techniques for mapping subsidence in oil fields, estimation of recharge to
shallow aquifers and freshwater lenses, calibration of satellite precipitation data and
mapping of hydrocarbon contamination using land surface temperature estimates.

Keywords EMR, hydrocarbon pollution, image enhancement, interferometry,

LST, TRMM

Contents
1 Fundamentals of Remote Sensing...................................................................................... 42
1.1 Electromagnetic Spectrum: The Photon and Radiometric Quantities ....................... 43
1.2 Electromagnetic Spectrum: Transmittance, Absorptance, and Reflectance.............. 45
2 Electromagnetic Spectrum: Distribution of Radiant Energies ........................................... 46
3 Electromagnetic Spectrum: Spectral Signatures ................................................................ 47
4 Sensor Technology ............................................................................................................. 48
4.1 Some Basic Image Processing Procedures ................................................................ 50
5 Processing and Classification of Remotely Sensed Data; Pattern Recognition;
Approaches to Data/Image Interpretation .......................................................................... 52
5.1 Digital Image ........................................................................................................... 52
5.2 Digital Image Analysis System ............................................................................... 53
5.3 Digital Data Formats ............................................................................................... 53

D. Al Ajmi and Saif ud din

Kuwait Institute for Scientific Research

© Springer-Verlag Berlin Heidelberg 2009

42 D. Al Ajmi and Saif ud din

5.4 Preprocessing .......................................................................................................... 54

5.5 Image Enhancement ................................................................................................ 55
5.6 Single Bands Linear Stretch .................................................................................... 57
5.7 Filtering ................................................................................................................... 58
5.8 Low-Frequency Filtering in the Spatial Domain ..................................................... 59
5.9 High-Frequency Filtering in the Spatial Domain .................................................... 61
5.10 Edge Enhancement in the Spatial Domain .............................................................. 61
5.11 Non-Linear Edge Enhancement .............................................................................. 63
5.12 Principal Components Analysis .............................................................................. 65
5.13 Band Ratio ............................................................................................................... 67
6 Classification ...................................................................................................................... 69
6.1 Image Classification ................................................................................................ 69
6.2 Unsupervised Classification .................................................................................... 70
6.3 Supervised Classification ........................................................................................ 74
6.4 Parallelepiped Classification ................................................................................... 75
6.5 The Minimum Distance to Mean Classification...................................................... 77
6.6 The Maximum Likelihood Classification................................................................ 77
7 Monitoring Applications of Environmental Consequences of War ................................... 79
7.1 Application I – Hydrocarbon Pollution consequent to 1991 Gulf War ................... 79
7.2 Application II – Precipitation Estimate Using Tropical Rainfall
Measuring Mission .................................................................................................. 90
7.3 Application III – Feature Extraction Technique for Palm Tree Census .................. 100
7.4 Application IV – Subsidence in Oil Fields.............................................................. 107
7.5 Application V – Sustainable Development of Fresh Water Resources ................... 113
References ................................................................................................................................ 120

1 Fundamentals of Remote Sensing

Remote sensing is a study of process, object, or phenomena without being in physi-

cal contact with it. The simplest example of remote sensing in daily life is viewing
the screen of your computer monitor (Fig. 1).
A physical quantity (light) emanates from the screen, which is a source of radia-
tion. The radiated light passes over a distance, and thus is “remote” to some extent,
until it is encountered and captured by the sensor (your eyes). Each eye sends a signal
to a processor (your brain) which records the data and interprets this into information.
Several of the human senses gather their awareness of the external world almost
entirely by perceiving a variety of signals, either emitted or reflected, actively or
passively, from objects that transmit this information in waves or pulses.
Remote sensing represents the acquisition and measurement of data/informa-
tion on some property(ies) of a phenomenon, object, or material by a recording
device not in physical, intimate contact with the feature(s) under surveillance;
techniques involve amassing knowledge pertinent to environments by measuring
force fields, electromagnetic radiation, or acoustic energy employing cameras,
radiometers and scanners, lasers, radio frequency receivers, radar systems, sonar,
thermal devices, seismographs, magnetometers, gravimeters, scintillometers, and
other instruments.
Remote Sensing 43

Fig. 1 Satellite image of Kuwait City

1.1 Electromagnetic Spectrum: The Photon

and Radiometric Quantities

The underlying basis for most remote sensing methods and systems is measuring
the varying energy levels of a photon. The variations in photon energies are tied to
the parameter wavelength or frequency. Electromagnetic radiation covers high to
low energy levels which combine into the electromagnetic spectrum (EMS).
Radiation from specific parts of the EMS contains photons of different wavelengths
whose energy levels fall within a discrete range of values. When any target material
is excited by internal processes or by interaction with incoming electromagnetic
radiation, it will emit photons of varying wavelengths whose radiometric quantities
differ at different wavelengths in a way diagnostic of the material.
The photon can be described as the messenger particle for EM force or as the
smallest bundle of light. This subatomic mass less particle comprises radiation
emitted by matter when it is excited thermally, or by nuclear processes (fusion, fis-
sion), or by bombardment with other radiation. It also can become involved as
reflected or absorbed radiation. Photons move at the speed of light: 299,792.46
kms−1. These particles also move as waves and hence, have a “dual” nature. These
waves follow a pattern that we described in terms of a sine (trigonometric) function,
as shown in two dimensions (Fig. 2).
44 D. Al Ajmi and Saif ud din

Fig. 2 Wave propagation and frequency in different bands

Wave length
Amplitude

de
Amplitu

Electric Wave

Magnetic Wave

Electromagnetic Spectrum: Transmittance, Absorptance, and Reflectance

Fig. 3 Electromagnetic spectrum: trandmittance, Absorptance, and Reflectance

The distance between two adjacent peaks on a wave is its wavelength. The total
number of peaks that pass by a reference in a second is frequency. A photon travels
as a wave having two components, oscillating as sine waves mutually at right
angles, one consisting of the varying electric field, the other the varying magnetic
field. Both have the same amplitudes and their max–min coincide. Unlike other
wave types which require a carrier (e.g., water waves), photon waves can transmit
through a vacuum (such as in space). When photons pass from one medium to
another, e.g., air to glass, refraction is observed (Fig. 3).
Remote Sensing 45

1.2 Electromagnetic Spectrum: Transmittance, Absorptance,

and Reflectance

Any beam of photons from some source passing through one medium to another
will experience any one or combination of any of the phenomena, i.e., transmission,
reflection, absorption, scattering.
The primary source of energy that illuminates natural targets is the Sun. The
main wavelength interval is between 2,000 and 34,000 Ångstrom [Å], with the
maximum power input close to 4,800 Å, which is in the visible green region. As
solar rays arrive at the Earth, the atmosphere absorbs or backscatters a fraction of
them and transmits the remainder. Upon striking the land and ocean surface, and
atmospheric targets, such as air, moisture, and clouds, the incoming radiation parti-
tions into three modes of energy-interaction response:
(1) Transmittance (t ) – some fraction of the radiation penetrates into certain sur-
face materials such as water and if the material is transparent and thin in one
dimension, normally passes through, generally with some diminution.
(2) Absorptance (a) – some radiation is absorbed through electron or molecular
reactions within the medium; a portion of this energy is then re-emitted, usually
at longer wavelengths, and some of it remains and heats the target;
(3) Reflectance (r) – some radiation reflects and scatters away from the target at
various angles, depending on the surface roughness and the angle of incidence
of the rays.
Because they involve ratios (to irradiance), these three parameters are dimen-
sionless numbers (between 0 and 1), but are commonly expressed as percentages.
Following the Law of Conservation of Energy: t + a + r = 1.
A fourth situation, when the emitted radiation results from internal atomic/molec-
ular excitation, usually related to the heat state of a body, is a thermal process.
There are two general types of reflecting surfaces that interact with EMR: specu-
lar (smooth) and diffuse (rough). These terms are defined geometrically, not physi-
cally. A surface may appear to be smooth in a physical sense, i.e., it appears and
feels smooth, but at a scale on the order of wavelengths of light, many irregularities
might occur throughout that surface. (A concrete roadway may appear smooth and
flat from a distance but feels rough when a finger passes over it, owing to small
grooves, pits, and protuberances.) Radiation impinging on a diffuse surface tends
to be reflected in many directions (scattered). The Rayleigh criterion is used to
determine surface roughness with respect to radiation.
A specular surface reflects radiation according to Snell’s Law which states that
the angle of incidence equals the angle of reflectance. Specular reflectances within
the visible wavelength range vary from as high as 0.99 for a very good mirror to as
low as 0.02–0.04 for a very smooth water surface.
In general, natural surfaces are almost always diffuse and depart significantly
from specular at shorter wavelengths (into the infrared) and may still be somewhat
diffuse in the microwave region.
46 D. Al Ajmi and Saif ud din

The term bidirectional reflectance describes the common observational condi-

tion in remote sensing in which the viewing angle j differs from the angle q of
incidence, and incoming/outgoing rays have different azimuths. Thus, reflectances
from the same target (type) change in value from various combinations of q and j:
this is particularly important when the sensor operates off-nadir and the Sun angle
and azimuth vary during the period of operation.

2 Electromagnetic Spectrum: Distribution of Radiant Energies

The electromagnetic radiation (EMR) extends over a wide range of wavelengths. A

narrow range of EMR extending from 0.4 to 0.7 mm, is detected by the human eye.
This is called the visible region. The EMR is divided into different regions based
on the wavelengths (Fig. 4).
At the very short wavelength end are gamma rays and X-rays normally meas-
ured in Ångstroms, which in the metric scale are in units of 10−10m. Radiation in
the ultraviolet extends from about 300 Å to about 4,000 Å. The visible region occu-
pies the range between 4,000 and 7,000 Å. The infrared region, spanning between
0.7 and 1,000 μm (or 1mm), has four subintervals of special interest: (1) reflected
IR (0.7–3.0 μm), and (2) its film responsive subset, the photographic IR (0.7–0.9
μm); (3) and (4) thermal bands at (3–5 μm) and (8–14 μm). The microwave region
spreads across 0.1–100 cm. The lowest frequency-longest wavelength region
beyond 100cm is the radio bands, from VHF (very high frequency) to ELF
(extremely low frequency); units applied to this region are often stated as frequencies
in units of Hertz (1Hz=1 cycle per second; KHz, MHz and GHz are kilo-, mega-,
and giga- Hertz, respectively).
The transmission of the energy bands through the atmosphere is variable, some
wavelengths are transmitted completely while others are absorbed completely.
Figure 5 depicts relative atmospheric radiation transmission and absorption at
different wavelengths.
Blue zones mark minimal passage of incoming and/or outgoing radiation,
whereas, white areas denote atmospheric transparency known as the window. Most
Gamma Rays

Microwave

Television
Ultraviolet

AM Ratio
Infrared
Visible
X-rays

Long Radio
Waves

10−14 10−12 10−10 10−8 10−6 10−4 10−2 1 102 104 106
Wave Length (m)
1 nm 1 um 1 cm 1 km

Fig. 4 Different regions of electromagnetic spectrum

Remote Sensing 47

V
I
S INFRARED MICROWAVE
I
B K X C S L P
L TRANSMISSION
14 13 12 11 10 RADAR BANDS 9 Frequency
E 10 10 10 FOR Q KM PATH 10 10 10 (Hz)
1.0 H2O
MOSTLY OPAQUE O2
o3 DUE TO H2O
CO2
H2O H2O O2
H2O CO2 H O O3 .5
H2O 2
CO2
H2O CO2 H2O H2O
H2O H2O

1.0 μm 10 100 μm 0.1 cm 1.0 cm 10 cm 100 cm

WAVE
LENGTH

Fig. 5 Atmospheric windows

remote-sensing instruments on board operate in one or more of these windows by

making their measurements with detectors tuned to specific wavelengths that pass
through the atmosphere. However, some sensors, especially those on meteorolo-
gical satellites, directly measure absorption phenomena, such as those associated
with carbon dioxide, CO2 and other gaseous molecules. The atmosphere is nearly
opaque to EM radiation in part of the mid-IR and all of the far-IR regions. In the
microwave region, by contrast, most of this radiation moves through unimpeded, so
radar waves reach the surface. Fortunately, absorption and other interactions occur
over many of the shorter wavelength regions, so that only a fraction of the incoming
radiation reaches the surface. Backscattering is a major phenomenon in the atmos-
phere. There are two prominent types of scattering Mie and Rayleigh. The Mie
scattering refers to reflection and refraction of radiation by atmospheric constituents
like smoke whose dimensions are comparable to the radiation wavelengths. Rayleigh
scattering results from constituents like atmospheric gases and water vapor that are
much smaller than the radiation wavelengths. Rayleigh scattering increases with
shorter wavelengths. The blue sky is a result of Rayleigh scattering.
Remote sensing of the Earth traditionally has used reflected energy in the visible
and infrared and emitted energy in the thermal infrared and microwave regions to
gather radiation that can be used to generate images whose variations represent dif-
ferent processes, phenomena, or objects. Images made from the varying wave-
length/intensity signals will show variations in gray tones in black and white
versions or colors (in terms of hue, saturation, and intensity) in colored versions.

3 Electromagnetic Spectrum: Spectral Signatures

The amount of solar radiation that reflects, absorbs, or transmits varies with wave-
length for any given material. This important property of matter makes it possible
to identify different substances or classes and separate them by their spectral signa-
ture. For example, at some wavelengths, sand reflects more energy than green
vegetation but at other wavelengths it absorbs more than does the vegetation. Using
reflectance differences, different crop types can be distinguished (Wheat; Alfalfa;
Potato; Tomato) (Fig. 6).
48 D. Al Ajmi and Saif ud din

0.6

0.4
Reflectence %

0.2 Reflectence curve

Wet sand
Dry sand
Wheat
Alfalfa
Potato
Tomato
0

0.4 0.6 0.8 1 1.2

Wave length (mm)

Fig. 6 Spectral signatures of different elements in variable wavelengths

4 Sensor Technology

Most remote sensing sensors are designed to measure frequency or wavelength of

the impinging radiation. The remote sensing sensors are grouped into several
classes according to their functionality and characteristics. A functional treatment
of several classes of sensors is plotted as a triangle diagram, in which the corner
members are determined by the principal parameter measured: Spectral; Spatial;
Intensity (Fig. 7).
Figure 8 presents a wider array of sensor types including active and passive sen-
sors of different types.
The Passive sensors acquire energy from an external source, e.g., the Sun, and
for Active the energy is generated from within the sensor system, beamed outward,
and the fraction returned is measured. Other types of commonly deployed sensors
in remote sensing studies are as follows:
Radiometer is an instrument that quantitatively measures the EM radiation in
some interval of the EM spectrum.
Spectrometer is an instrument with a component, such as a prism or diffraction
grating, which can break radiation extending over a part of the spectrum into dis-
crete wavelengths and/or separate them.
The term spectroradiometer tends to imply that the dispersed radiation is in
bands rather than discrete wavelengths. Most of the remote sensing satellites have
spectroradiometers on board.
Remote Sensing 49

Fig. 7 Different devices based on principal parameters they measure

There are many remote sensors

Microwave radiometer
non imaging
Magnetic sensor
non-scanning Monochrome
Gravimeter
Natural Color
Fourier spectrometer
Passive Infrared
Other (Resistivity, etc.)
imaging Color Infrared
Camera
Others
Sensor TV camera
type Image plane scanning
scanning imaging Solid scanner
Optical mechanical
Object plane scanning scanner
Microwave radiometer Microwave radiometer
non-scanning non-imaging
Microwave altimeter
Active Laser water depth meter
Laser distance meter
Image plane scanning Passive phased array radar
scanning imaging
Real aperture radar
Object plane scanning
Synthetic aperture radar

Fig. 8 Types of Sensors

50 D. Al Ajmi and Saif ud din

Sensors that instantaneously measure radiation coming from the entire scene at
once are called framing systems, i.e., eye, camera, and TV. The size of the scene
that is framed is determined by the apertures and optics in the system defines the
field of view.
If the scene is sensed point by point along successive lines over a finite time, this
mode of measurement is a scanning system. Most non-camera sensors operating
from moving platforms image the scene by scanning.
The radiation – normally visible and/or Near and Short Wave IR, and/or thermal
emissive in nature – must then be broken into its spectral elements, into broad to
narrow bands. The width in wavelength units of a band or channel is defined by the
instrument’s spectral resolution. A more vital aspect of sensor characteristic and
performance is spatial resolution. Spatial resolution represents the “ability to rec-
ognize and separate features of specific sizes.” The common definition of spatial
resolution is often simply stated as the smallest size of an object that can be picked
out from its surrounding objects or features. This separation from neighbors or
background may or may not be sufficient to identify the object.
Three variables control the achieved spatial resolution: (1) the nature of the
target features, the most important being size; (2) the distance between the target
and the sensing device; and (3) some inherent properties of the sensor embodied in
the term resolving power.
The spectral properties of the material help in differentiating them from one
another. Typical spectra are shown in Fig. 9 the figures show reflectance for veg-
etation rises abruptly at about 0.7 μm, followed by a gradual drop at about 1.1 μm.
The first spectral signatures indicate a gradual rise in reflectance with increasing
wavelengths for those particular common manmade materials on the ground.
Concrete, being light-colored and bright, has a notably higher average than dark
asphalt. The other materials fall in between. The shingles are probably bluish in
color as suggested by a rise in reflectance from about 0.4–0.5 μm and a flat
response in the remainder of the visible (0.4–0.7 μm) light region. The second
curves (on the right or bottom) indicate most vegetation types are very similar in
response between 0.3 and 0.5 μm; show moderate variations in the 0.5–0.6 μm
interval; and display maximum variability (hence optimum discrimination) in the
0.7–0.9 μm range (Fig. 9).

4.1 Some Basic Image Processing Procedures

The basic image processing involves computer-based processing procedures in

highlighting and extracting information about scene content, that is, the recogni-
tion, appearance, and identification of materials, objects, features, and classes. The
processing procedures can broadly be characterized into three categories: image
restoration (preprocessing); image enhancement; and classification and informa-
tion extraction.
Remote Sensing 51

a NONVEGETATED LAND AREAS

70 CONCRETE
ASPHALT
FROM ROOT AND MILLER (1971)
BARE SOIL
60
GRAVEL
PERCENT REFLECTANCE

SHINGLES
50

40

30

20

10

0
.3 .4 .5 .6 .7 .8 .9 1.0 1.1 1.2 1.3
WAVELENGTH (μm)

b VEGETATED LAND AREAS

70 GRASS
TREE
60 SUGAR BEET
WHEAT STUBBLE
PERCENT REFLECTANCE

FALLOW FIELDS
50

40

30

20

10

0
.3 .4 .5 .6 .7 .8 .9 1.0 1.1 1.2 1.3
WAVELENGTH (μm)

Fig. 9 Spectral signatures of vegetated and non vegetated areas in different wavelengths
52 D. Al Ajmi and Saif ud din

5 Processing and Classification of Remotely Sensed Data;

Pattern Recognition; Approaches to Data/Image
Interpretation

The remotely sensed data from the satellites were analyzed for identifying objects,
extraction of information on surface features and deducing their properties through
observations made on the reflected/scattered energy from the earth in different
spectral bands. The images recorded in digital form are processed on the computer
to produce images for interpretation. The digital images contain large amounts of
data and computers with large data handling capabilities are used for digital image
processing (DIP) work.
Digital image processing involves techniques for manipulation of digital image
data by computers to generate analog signals. Digital image processing comprises
the operations for noise removal, geometric and radiometric corrections, enhancement
of images, data compaction, image display and recording, image data manipulation
and management.
In India NRSA provides data, which are corrected for noise, geometric and
radiometric distortions. The digital image processing operations carried out during
the present study involved conversion of data from band interleaved by line (BIL)
to intergraph format (.cot) and numerical operations for image enhancement, infor-
mation extraction, data compaction, image display and recording and image data
manipulation and management.

5.1 Digital Image

Digital image consists of discrete picture elements called pixels. Associated with
each pixel is a number represented as Dn (Digital number), that depicts the average
radiance of a relatively small area within a scene (pixel). The size of this area deter-
mines the reproduction of details within the scene. As the pixel size is reduced,
more scene details are preserved in a digital representation. The pixel has a gray
scale value where 0 corresponds to black and 255 for white in between there are
shades of gray [1–5].
A Dn is simply one of a set of numbers based on powers of 2, such as 26 or 64.
The range of radiances, which instrument-wise, can be, for example, recorded as
varying voltages if the sensor signal is one which is, say, the conversion of photons
counted at a specific wavelength or wavelength intervals. The lower and upper limits
of the sensor’s response capability form the end members of the Dn range selected.
The voltages are divided into equal whole number units based on the digitizing range
selected. Thus, a Landsat TM band can have its voltage values – the maximum and
minimum that can be measured – subdivided into 28 or 256 equal units. These are
arbitrarily set at 0 for the lowest value, so the range is then 0–255.
Remote Sensing 53

5.2 Digital Image Analysis System

The digital image processing system comprises hardware and software elements
which help the analyst in extracting meaningful information from the data. The basic
system components include input devices, processing elements, interactive devices
and output devices. The DIP systems are capable of processing large volumes of data
at a very fast rate. The speed of operation is of immense importance in DIP as some
of the image analyses techniques like classification, spatial filtering, etc. are highly
computation bound, to facilitate this fast processing high speed pipe line processors
like array processors may be used in conjunction with the CPU.
A display system is an interactive device which is necessary for the user to inter-
act with the computer. Normally a color scrolling display system with an ability to
zoom, shrink and pan is provided. The display area could either be 512×512
(pixels×scan lines), 1,024×1,024 or even more in case of high-resolution display
systems. The state-of-the-art display systems are microprocessor-controlled and
provide a wide range of image analysis capabilities at an increased speed. The
analyzed outputs are to be stored in a form which aids further analysis by the users.
Outputs can be any of the following:
• Computer compatible tapes (CCTs)/floppies/cartridges/CDs
• Gray scale maps on line printer or color plotter outputs
• Photographs of color monitor displays
• High precision film output/imagery

5.3 Digital Data Formats

The remotely sensed data acquired from the satellites are stored in different types
of formats.
• Band sequential (BSQ)
• Band interleaved by line (BIL)
• Band interleaved by pixel (BIP).
Each of these formats is preceded on tape by header or tailer information, which
consists of ancillary data about the date, altitude of the sensor, sun angle, and so on.
This information is useful when the data is corrected geometrically or radiometri-
cally. The data are normally recorded on nine-track CCTs with data density on the
tape of 800, 1,600, or 6,250 bits per inch (bpi).

5.3.1 Band Sequential Format (BSQ)

This format requires that all data for a single band covering the entire scene be writ-
ten as one file. Thus, if one wants the area in the center of the scene in four bands,
54 D. Al Ajmi and Saif ud din

it would be necessary to read into this location in four separate files to extract the
desired information. The format is preferred as it is not necessary to read serially
past unwanted information if certain bands are of no value. The number of tapes may
be dependent on the number of bands provided for the scene.

5.3.2 Band Interleaved by Line Format (BIL)

In this format the data for bands are written line by line onto the tape (i.e., line 1 band
1, line 1 band 2, line 1 band 3, line 1 band 4). It is a useful format if all the bands are
to be used in the analysis. If some bands are not of interest, this format will be inef-
ficient, since it is necessary in BIL to read serially past all the unwanted data.

5.3.3 Band Interleaved by Pixel Format (BIP)

In this format, the data for the pixel in all bands are written together. Taking the
example of LANDSAT – MSS 4 band data every element in the matrix has 4 pixel
values pixel 1,1, of band 1; pixel 1,1, of band 2; pixel 1,1 of band 3; pixel 1,1 of
band 4, pixel 1,2 of band 1, pixel 1,2 of band 2, pixel 1,2 of band 3; pixel 1,2 of
band 4; etc. This data format is of use if all the bands are to be used, otherwise it
will be inefficient. This format is not very popular now; it was extensively used by
the EROS Data Center for LANDSAT scenes in initial stages.

5.4 Preprocessing

Preprocessing is an important and diverse set of image preparation programs that

act to offset problems with the band data and recalculate Dn values that minimize
these problems. Among the programs that optimize these values are atmospheric
correction; sun illumination geometry; surface-induced geometric distortions;
spacecraft velocity and attitude variations; effects of Earth rotation, elevation, cur-
vature, abnormalities of instrument performance; loss of specific scan lines, and
others. Once performed on the raw data, these adjustments require appropriate
radiometric and geometric corrections.
Resampling is one approach commonly used to produce better estimates of the
Dn values for individual pixels. After the various geometric corrections and transla-
tions have been applied, the net effect is that the resulting redistribution of pixels
involves their spatial displacements to new, more accurate relative positions.
However, the radiometric values of the displaced pixels no longer represent the real
world values that would be obtained if this new pixel array could be resensed by
the scanner. The particular mixture of surface objects or materials in the original
pixel has changed somewhat. Three most common transformations are:
Remote Sensing 55

The nearest neighbor technique, the transformed pixel takes the value of the
closest pixel in the pre-shifted array.
The bilinear interpolation approach, the average of the Dns for the 4 pixels
surrounding the transformed output pixel is used.
The cubic convolution technique averages the 16 closest input pixels; this
usually leads to the sharpest image.

5.5 Image Enhancement

The principal objective of the enhancement techniques is to process a given

image so that the result is better than the original image. Image enhancement
techniques improve the quality of the image. These techniques are most useful as
many satellite imageries when examined on a color display give inadequate infor-
mation for image interpretation. There exists a wide variety of techniques for
improving image quality. Normally image enhancement involves techniques for
increasing the visual distinction between features in a scene. Basically new
images are created from the original image data in order to increase the interpret-
ability of information that can be visually extracted from the data, the following
techniques were used:
• Contrast enhancement
• Density slicing
• Edge enhancement
• Spatial filtering.

5.5.1 Contrast Enhancement

The sensors record reflected and radiant flux from the earth’s surface materials. The
reflectance of materials in different bands of electromagnetic radiation varies, this
results in contrast between the materials when recorded by remote sensing systems.
Sometimes different materials reflect similar amounts of radiant flux in certain
bands of electromagnetic spectrum resulting in a relatively low contrast image.
Besides low contrast characteristic of materials, the lowering of the sensitivity of
detectors often results in low contrast imagery.
The contrast can be defined as the ratio of the maximum intensity to the mini-
mum intensity over an image.

C = Imax/ Imin.

Contrast ratio has a strong bearing on the resolving power and detectability of
an image. The larger this ratio, the easier it is to interpret the image outline. Most
of the satellite images lack adequate contrast and require contrast improvement.
Low contrast may result from the following:
56 D. Al Ajmi and Saif ud din

• The object and the background of the terrain have nearly uniform electromag-
netic response in the wavelength band of energy that is recorded by the remote
sensing system, i.e., the scene itself has a low contrast.
• Scattering of electromagnetic energy by the atmosphere can reduce the contrast
of a scene. This effect is most pronounced in the shorter wavelength portions of
EMR.
• The remote sensing system may lack sufficient sensitivity to detect and record
the contrast of the terrain. Also, incorrect recording techniques can result in low
contrast imagery although the scene has a high contrast.
The imageries with low contrast are commonly referred to as “washed out” with
nearly uniform tones of gray. Detectors on board are designed to record a wide
range of scene brightness values without getting saturated. The image of the study
area has not utilized the full range of brightness, unless it is subjected to
enhancement.
The contrast enhancement technique expands the range of brightness values in an
image so that the image can be efficiently displayed in a desired manner. The Dn
values are literally pulled apart, that is, expanded over a full range of 0–255 Dn val-
ues. The effect is to increase the visual contrast between two areas of initially nearly
uniform Dn values, which provide ease in identification. The contrast modification is
the most commonly applied image enhancement technique. Some enhancement is
possible by photographic duplication on a high contrast film. This photographic
enhancement results in an overall loss of information. The limited dynamic range of
film results in loss of information at the bright and dark extremes of image to avoid
loss of information at the tails, the gray scale is enhanced by digital image
enhancement.
The contrast enhancements are of two types, i.e., linear and non-linear contrast
enhancement.

5.5.2 Linear Contrast Enhancement

The linear contrast enhancement expands the original Dn values to make use of the
total range of 256 gray levels of the output device. Linear contrast enhancement is
best applied to remotely sensed images with Gaussian or near Gaussian histograms,
that is when all the brightness values fall between a single narrow range and only
one mode is apparent [2, 5]. To perform linear contrast enhancement the image
statistics are examined to determine minimum and maximum brightness values in
the band (mink and maxk). The output brightness value (BVout), is computed according
to the following equation:

BVout = (BVin − min k ) / (min k − max k ) quant k ,

where BVin is the original input brightness value and quantk is the range of bright-
ness value that can be displayed on CRT (e.g., 256).
Remote Sensing 57

5.5.3 Min–Max Contrast Stretch

The minimum and maximum Dn values are redistributed in this contrast enhance-
ment from 0 to 255 keeping the min at 0 and max at 255.

5.5.4 Percentage Linear Contrast Stretch

The mink and maxk that lie in a certain percentage of pixels can be specified from
the mean of the histogram. This is known as percentage linear contrast stretch. In
this enhancement if (minimum–maximum) the standard deviation is ±1 and mini-
mum and maximum Dn values are 20 and 60, then on stretching the Dn value 20
will become 0 and 60 will be at 255 and all the values from 0 to 19 will be 0 and
all the values from 61 to 255 will be 255. This results in more pure black and white
pixels in the scene, increasing the contrast of the image but the information content
of pixels that are saturated at 0 and 255 is lost. The slope of a percentage linear
contrast stretch is much greater than for a simple min–max stretch [2, 5].
The linear stretch of an image having a non-Gaussian histogram, is performed
piece wise, here a number of linear segments are identified and are stretched.

5.6 Single Bands Linear Stretch

5.6.1 Non-Linear Contrast Enhancement

In this type of enhancement, the input and output data values follow a non-linear
transformation.

5.6.2 Gaussian Stretch

The brightness levels of the individual bands are redistributed to an approximate

normal distribution. The near approximate normal distribution of the histogram
stretching is achieved by changing mean and standard deviation of the stretched
histogram. The mean of the Dn values of the histogram is adjusted by shifting the
center of the normal distribution in the range of 0–255 Dn. The standard deviation
varies between 1 and 255 Dn. The increase in standard deviation widens the normal
distribution curve; it will increase the contrast of the image. The decrease of stand-
ard deviation results in a bunch up of stretched histogram about the mean, which
reduces the contrast of the image.

5.6.3 Histogram Equalization

One of the most useful techniques is histogram equalization. In this the algorithm
passes through the individual bands of the dataset and assigns approximately an
58 D. Al Ajmi and Saif ud din

equal number of pixels to each of the user-specified gray scale. The histogram
equalization provides greatest contrast enhancement to the most populated range of
brightness values in the image. It automatically reduces contrast in the very light or
dark parts of the image associated with tails of a normally distributed histogram.
The frequency of occurrence of brightness values, f(BVi) will be the ratio of the
number of pixels in the scene with the same value to the total number of pixels in
the scene. The probability pi=f(BVi)/n, where p is probability and n is total number
of pixels. For each value level BVi in the quant range of 0–255 of the original
histogram, a new cumulative frequency Ki is calculated:
quant k

Ki = ∑
i=0
f (BVi ) / n,

where the summation counts the frequency of pixels in the image with brightness
values equal to or less than BVi and n is total number of pixels in the entire scene.
The histogram transformation function Ki with the original values Li to determine
the closest match and is reassigned to an appropriate brightness value. The histo-
gram equalization results in rescaling of brightness levels of the scene into a lower
number of brightness levels. The quantization of the gray levels reduces the
smoothness of the edges.

5.6.4 Logarithmic Contrast Enhancement

The logarithmic contrast enhancement is non-linear; this enhancement has the

greatest impact on the brightness values found in darker parts of the histogram. It
could be reversed to enhance values in the brighter part of the histogram by scaling
the input data using an inverse log function.

5.7 Filtering

The filtering of the data can be carried out both in spatial and frequency domains.
The spatial filters are used as masks. The enhancement in the frequency domain
involved the Fourier Transform of the image and its multiplication by a filter trans-
fer function. The inverse transform of the result can be taken to produce enhanced
images in the frequency domain. The most common spatial filters include lowpass,
highpass, and median filters.

5.7.1 Spatial Filtering

Spatial Frequency is defined as the number of changes in the brightness value per
unit distance or any particular part of the image. If there are very few changes in
brightness value over a given area in an image, it is called a low frequency image
Remote Sensing 59

and if the changes are frequent over a short distance it is called a high frequency
image. Spatial frequency describes the brightness values by taking into account the
brightness of the neighboring pixels. The spatial frequency in an image can be
spatially enhanced or subdued using two different approaches, i.e., spatial convolution
(spatial filtering) and Fourier analysis (frequency domain). The spatial convolution
is used for enhancing the low and high frequency detail in imagery.

5.8 Low-Frequency Filtering in the Spatial Domain

5.8.1 Low-Pass Filters

These are enhancements that block the high spatial frequency detail. They are
called low-pass filters. This filter evaluates a particular input pixel brightness value
BVin, and the pixel surrounding the input pixel and output is a new brightness
value (BVout) for the central pixel, that is the mean of this convolution. The size of
the kernel (n) is usually 3×3, 5×5 and 7×7.
If the coefficients in a low frequency mask are set equal to 1

1 1 1
Mask A = 1 1 1
1 1 1

The coefficients ci, in the mask are multiplied by the following individual bright-
ness values (BVi) in the input image

c1 × BV1 c2 × BV2 c3 × BV3

Mask template = c4 × BV4 c5 × BV5 c6 × BV6 ,
c7 × BV7 c8 × BV8 c9 × BV9

where BV1 = BVi−1,j−1 BV2 = BVi−1,j BV3 = BVi−1,j+1 BV4 = BVi,j−1 BV5 = BVi,j BV6=BVi,j+1
BV7 = BVi+1,j−1 BV8 = BVi+1,j BV9 = BVi+1,j+1.
The primary pixel under investigation is BV5=BVi,j. The original data which will
result in a low frequency image is given by expression

n
Low - Frequency Image = Int.∑ ci × BVi / n
i =1

= Int.(BV1 + BV2 + BV3 +  + BV9) / 9.

The values for each pixel are calculated and this is called image smoothing. The
image smoothing is useful in removing periodic noise recorded by remote sensing
systems. The effect of low-pass filtering is smoothing the image by cutting of the
high frequency components and allowing only low frequency.
60 D. Al Ajmi and Saif ud din

5.8.2 Weighted Low-Pass Filters

The weighted low-pass mask carries an option of giving a center weight to an indi-
vidual pixel. It can be used to have selective smoothing.

5.8.3 Median Filter

The median filter is useful for removing noise in an image, especially shot noise
by which individual pixels are corrupted or missing, instead of taking the mean of
nine pixels in a 3×3 mask, the median filter ranks the pixels in the neighborhood
from lowest to highest and selects the median value. The median filter has several
advantages when compared with weighted convolution filters. It does not shift
boundaries, has minimal degradation to edges, allows the median filter to be applied
repeatedly, as a result fine details are erased and large regions acquire the same
brightness value. The standard median filter will erase some lines in the image
which are narrower than the half-width of the neighborhood and round or clip
corners.

5.8.4 Edge-Preserving Median Filter

The edge preserving median filter is where the median value of the black pixels and
gray pixels is computed in a 5×5 array, these two values and the central original
brightness value are ranked in ascending order. A final median value is selected to
replace the central pixel. This filter preserves edges and corners.
The minimum-maximum filter operates on one pixel at a time and examines
the brightness value in a user-specified radius (3×3 kernel) and replaces bright-
ness values of the current pixel with the minimum or maximum brightness value
encountered.

5.8.5 Sigma Filter

An adaptive box filter (Sigma) is of value for removing noise in digital images. The
adaptive box filters are used to remove random bit errors where pixel values have
no relation to image scene, i.e., shot noise and for smooth noisy data where the
pixels are related to image scene but with an additive or multiplicative component
of noise.
The procedures rely on computation of the standard deviation components s of
only those pixels within a local box surrounding the central pixel [i.e., the eight
values surrounding the fifth (central) pixel in a 3×3 mask]. The original brightness
value at location 5 is considered to be a bit error if it deviates from the box mean
of eight values by more than 1.0–2.0s and is replaced by the box mean. This is
called the adaptive filter because it is based on computation of standard deviation
Remote Sensing 61

for each 3×3 window rather than on standard deviation of the entire scene. In this
very minor bit errors are removed from low variance areas, but valid data along
sharp edges and corners are not replaced.

5.8.6 Lee Filter

The Lee adaptive filter can be used for cleaning up of extremely noisy images and
is based on the sigma filter [6]. Lee’s filter computes standard deviation for the
entire scene. Then the fifth pixel in a 3×3 moving window is replaced by the aver-
age of the neighboring pixels that had the intensity within a fixed s range of the
central pixel. The filter averaged only those pixels within the box that had intensi-
ties within 1.0–2.0s of the central pixel. This technique effectively reduces speckle
and the salt and pepper texture in images without eliminating fine details [7]. The
two filters can be combined into a single program for processing images with both
random bit error and noisy data.

5.9 High-Frequency Filtering in the Spatial Domain

High pass filtering is applied to imagery to remove the slowly varying components
and enhance the high frequency local variations. High frequency filter (HFF) is
computed by subtracting the output of the low frequency filter (LFF) from twice the
value of the original central pixel in a matrix of 3×3, 5×5, or 7×7. In a 3×3 matrix
the output value of the fifth pixel is computed by the following expression:

HFF5 out = (2 × BV5 ) − LFF5 out .

Brightness values tend to be highly correlated in a 3×3 window, thus the highest
frequency filtered images will have a relatively narrow intensity histogram. Thus,
the output from most high-frequency filtered images must be contrast stretched
prior to visual analysis. This high pass filter sharpens edges.

5.10 Edge Enhancement in the Spatial Domain

Edge enhancement delineates the edges and makes the shapes and details compris-
ing the image more conspicuous and perhaps easier to analyze. Generally, the picto-
rial edges are simply sharp changes in brightness value between two adjacent
pixels. In earth science applications the most valuable information that may be
derived from an image is contained in the edges surrounding various objects of
interest. The edges may be enhanced using either linear or non-linear edge-
enhancement techniques.
62 D. Al Ajmi and Saif ud din

5.10.1 Linear Edge Enhancement

A straightforward method of extracting edges in remotely sensed imagery is appli-

cation of a directional first-difference algorithm that approximates the first derivative
between two adjacent pixels. The algorithm produces the first difference of
the image input in the horizontal, vertical, and diagonal directions. The algorithm
for enhancement is:
Vertical: BVi,j=BVi,j−BVi,j+1+K
Horizontal: BVi,j=BVi,j−BVi–1,j+K
NE diagonal: BVi,j=BVi,j−BVi+1,j+1+K
SE diagonal: BVi,j=BVi,j−BVi–1,j+1+K
The result of the subtraction can either be negative or positive, a constant K,
usually 127 is added to make all values positive and centered between 0 and 255.
This enhances contrast of adjacent pixels. The resultant image is normally min–
max contrast stretched to enhance the edges even more. It is best to make the min–
max values in the contrast stretch a uniform distance from midrange value (127).
This causes the uniform areas to appear in shades of gray, while the important edges
become black or white. The edge enhancement is carried out by convoluting the
original data with a kernel. Chavez et al. [8] suggested the optimum kernel size as
3×3, 5×5, or 7×7 based on the surface roughness and sun angle characteristics of
the data. They developed a procedure based on the first difference in the horizontal
direction. An offset of 127 is added to the result and the data contrast stretched. The
direction of the embossing is controlled by changing the coefficient around the
periphery of the mask.

5.10.2 Gradient Filter

It had been observed that low-pass filtering based on averaging of pixel values
resulted in smoothing; the process is analogous to integration. The differentiation
can be expected to have an opposite effect, i.e., sharpening of the image. The object
is not recognized in the image because of the difference of gray values but also
because of the difference in the pattern and orientation of pixels. The gradient filter
was developed based on the logic that the local orientation of a pattern is the
property which describes the edges of the image features.
In an image for a function ƒ(x,y), the gradient of ƒ at coordinate (x,y) can be
defined as vector
⎡ df dx ⎤
∇F = ⎢ ⎥.
⎣ df dy ⎦
The magnitude of this vector,

∇f = mag(∇f ) = [(df / dx )2 + (df / dy)2 ]1/ 2 .

Remote Sensing 63

5.10.3 Compass Filter

The compass mask is based on gradient filtering which may be used to perform
two-dimensional, discrete differential directional edge enhancement. The compass
filter suggests the slope direction of maximum response. Thus, the east gradient
mask produces a maximum output for horizontal changes in brightness value from
west to east. The gradient mask having zero weight, results in no output response
over regions with constant brightness values, i.e., where no edges are present.

5.10.4 Laplacian Filter

The Laplacian filter is a special high pass filter, effectively enhancing the plume and
other subtle sensor noise in the image. A Laplacian convolution mask is applied to
perform edge enhancement and is insensitive to direction and invariant to rotation.

5.10.5 The Laplacian Filter States

Dx 2 = − Dx + Dx ,
L = Dx 2 + Dy 2 ,
 ) + 2 Cos (pKy
Lˆ = 2 Cos (pKx  ) − 4.

The subtraction of the Laplacian edge enhancement from the original image
restores the overall gray variation. It also sharpens the image by locally increasing
the contrast at discontinuities. The Laplacian operator highlights points, lines, and
edges in the image and suppresses uniform and smoothly varying regions. By itself,
the Laplacian image is difficult to interpret.
The combination of gradient and Laplacian edge operators may be used for edge
enhancement, which may be superior to either edge enhancement alone.

5.11 Non-Linear Edge Enhancement

Non-linear edge enhancements are performed using non-linear combinations of

pixels. Many algorithms are applied using either 2×2 or 3×3 kernels.

5.11.1 Sobel Edge Detector

The Sobel edge detector is based on the notion of the 3×3 window and is computed
according to the relationship:

Sobel 5, out = √ x 2 + y 2 .
64 D. Al Ajmi and Saif ud din

Here X = (BV3 + 2BV6 + BV9 )–(BV1 + 2BV4 + BV7), Y = (BV1 + 2BV2 + BV3)–
(BV7+2BV8+BV9 ).
The operator detects horizontal, vertical, and diagonal edges. Each pixel in an
image is declared an edge if its Sobel values exceeds some user-specified threshold.
This information may be used to create the edge map, which often appears as a
white line on a black background or vice versa.

5.11.2 Robert Edge Detector

Robert’s edge detector is based on the use of only four elements of 3×3 mask [9].
The new pixel value at pixel location BV5,out is computed according to equation:

Roberts5, out = x + y,

where X= |BV5–BV9|, Y= |BV6–BV8|.

5.11.3 Kirsch Filter

The gradient filters used to imply enhancement of edge gradient in two orthogonal
directions, i.e., row and columns. In Kirsch the edge enhancement is carried out in
8 modulo. The gradient directions are East, North-east, North, North-west, West,
South-west, South, and South-east. The template gradient 3×3 impulse response
arrays were used where the scale factor was 1/15 [10].
The Kirsch non-linear edge enhancement calculates the gradient at pixel loca-
tion BVi,j. The algorithm applied is:

{ 7

}
BVi,j = max 1,max [Abs(5Si – 3Ti)] ,
i=0

where Si=BVi+BVi+1+BVi+2 and Ti=BVi+3+BVi+4+BVi+5+BVi+6+BVi+7.

The subscripts of BV are evaluated modulo 8, meaning that the computation
moves around the perimeter of the mask in eight steps. The edge enhancement
computes the maximal compass gradient magnitude about input image point BVi,j.
The value of Si equals the sum of three adjacent pixels, while Ti equals the sum of
the remaining four adjacent pixels. The input pixel value at BVi,j is never used in
computation. The effect of the Kirsch operator is enhancement of linear edges and
boundaries.

5.11.4 Wallis Filter

The edge crisping has also been achieved through statistical differentiation which
involves generation of an image by dividing each pixel value by its estimated stand-
ard deviation according to the basic relation:
Remote Sensing 65

G( j , k ) = F( j , k ) S( j , k ) ,

where F(j,k) of the pixel S(j,k) is the expected standard deviation and G(j,k) is the com-
puted Dn value. The statistical differentiation operator in which the enhanced
image takes values from desired first-order and second-order moments is referred
to as the Wallis filter [11].

5.11.5 Prewitt Filter

The Prewitt filter is an edge gradient operator described by pixel numbering con-
vention [12]. In the Prewitt filter the row and column gradients are normalized to
provide unit gain and positive weighted and unit gain negative weighted averages
about a separated edge position. Unlike the Sobel edge detector the Prewitt edge
detector has K=1, whereas the Sobel detector has K=2, and as a result the pixel
values of the north, south, east, and west remain the same. The output of the
Prewitt Filtering usually eliminates the edges in the scene as compared to Sobel
output, where the north, south, east, and west pixel values are doubled.

5.11.6 Frost Filter

The frost command filter removes high frequency noise while preserving the edges.
The command removes “speckle” noise from images, besides it can be used to
remove high frequency noise from any type of image. The areas of the images having
low spatial frequency are smoothed and the areas containing high spatial frequen-
cies are not affected. The result is the speckle noise is removed from smooth areas
of the image while the edges are kept clean.

5.12 Principal Components Analysis

The principal component analysis (PCA) or Karhunen–Loeve analysis is useful for

analysis of highly correlated multi-spectral remotely sensed data [13, 14]. The
transformation of raw remote sensor data using PCA can result in new principal
component images that may be more interpretable than the original data [14, 15].
For PCA the transformation is applied to a correlated set of multi-spectral data,
application of the transformation to the correlated remote sensor data will result in
another uncorrelated multi-spectral dataset that has certain ordered variance proper-
ties. This transformation is conceptualized by considering the two-dimensional
distribution of pixel values obtained in two bands that can be labeled as X1 and X2.
The spread or variance of the distribution of points is an indication of the correla-
tion and quality of information associated with both bands, if all the points are
clustered in an extremely tight zone in two-dimensional space, these data will
provide very little information.
66 D. Al Ajmi and Saif ud din

The initial measurement coordinate axes X1 and X2 may not be the best arrange-
ment in multi-spectral feature space to analyze the remotely sensed data associated
with these two bands. The PCA will translate and/or rotate the original axes so that
the original brightness values on axes X1 and X2 are redistributed (re-projected) onto
a new set of axes as X’1 and X’2. The X’1 coordinate system might then be rotated
about its new origin (m1,m2) in the new coordinate system some f degrees so that the
first axis X’1 is associated with maximum amount of variance in the scatter point,
this new axis is called the first principal component axis PC1=l1. The second prin-
cipal component axis PC2=l2 is orthogonal to PC1. Thus, the major and minor axis
of ellipsoid of points in bands X1 and X2 are called the principal components. To
re-project the original data on the X1 and X2 axes onto the PC1 and PC2 axes certain
linear transformations are to be applied to the original pixel values. The linear
transformation required is derived from the covariance matrix of the original data
set. Thus, this is a data-dependent process with each data set yielding different
transformation coefficients. The transformation is calculated from the original
spectral statistics [16] as follows:
• The n×n covariance matrix, covariance of the n-dimensional remote sensing data
set to be transformed is computed. Use of the covariance matrix results in an
unstandardized PCA, whereas use of the correlation matrix results in a standard-
ized PCA.
• The eigenvalues E=[l1,1, l2,2, l3,3,…….., ln,n] and eigenvectors EV=[akp… for
k=1 to n bands and p=1 to n components].

EVT = EV Cov,
where EVT is the transpose of the eigenvector matrix, EV, and E is a diagonal
covariance matrix whose elements lii, called eigenvalues, are the variance of the
pth principal components, where p=1 to n components. The non-diagonal eigen-
values, lii are equal to zero and therefore can be ignored. The number of non-
zero eigen values in an n×n covariance matrix always equal n, the number of
bands examined. The eigen values are often called components (eigen value 1
may be referred to as PC1) where lp is the pth eigenvalue out of the possible n
eigen values. By calculating the correlation of each band k with each component
p, it is possible to determine how each band is associated with each principal
component. The equation is

Rkp = akp × √ lp √ Vark ,

where akp, eigenvector for band k and component p; lp, pth eigenvalue; Vark=variance
of band k in the covariance matrix.
The above computation results in a new n×n matrix. Each component contrib-
utes different information and it can be seen in different images. To do this it is
necessary to identify the brightness values (BVi,j,k) associated with each pixel.
Loeve [17] and Castleman [18] suggest that standardized PCA based on compu-
tation of eigenvalues from correlation matrices is superior to unstandardized PCA
Remote Sensing 67

computed from covariance matrices when analyzing change in multi-temporal

image datasets.
Another use of PCA, mainly as a means to improve image enhancement, is
known as a decorrelation stretch (DS). The DS optimizes the assignment of colors
that bring out subtle differences not readily distinguished in natural and false color
composites. This reduction in interband correlations emphasizes small but often
diagnostic reflectance or emittance variations owing to topography and tempera-
ture. The first step is to transform band data into at least the first three PCs. Each
component is rescaled by normalizing the variance of the PC vectors. Then each PC
image is stretched, usually following the Gaussian mode. The stretched PC data are
then projected back into the original channels which are enhanced to maximize
spectral sensitivity (Fig. 10).
Users of ASTER data have found decorrelation stretching to be particularly
effective in image display. The stretch is effective whether the bands used are in the
Visible, the SWIR, or the thermal IR interval. The three ASTER scenes in Fig. 11
(again, of an unidentified area) show the effects of a DS.
The difference between the PC color composites and the DS color composites is
generally not large, but extra statistic data manipulation in the latter often leads to
a better product.

5.13 Band Ratio

The band ratio minimizes the effect of environmental factors like slope and aspect,
shadows or seasonal changes which affect the brightness values of an identical
surface. The band ratios provide unique information that is not available in any
single band. The mathematical expression of the ratio function is:

BVi , j ,r = BVi , j , k BVi , j ,l ,

Fig. 10 An Original image (right) and PC enhanced imagery (left)

68 D. Al Ajmi and Saif ud din

Fig. 11 Decorrelation stretch of visible, SWIR and thermal IR images from ASTER

where BVi,j,r is the output ratio value for a pixel at row i, column j, ratio of band r;
BVi,j,k is the brightness value at the same location in band k and BVi,j,l is the brightness
value in band l. If the value of brightness BVi,j is equal to 0, then the alternatives
are to be taken as, the mathematical domain of the function as 1/255 to 255 through
0; or to assign 0 a value of 1; or a small value of the order of 0.1 can be added
to the denominator.
Normalization is done to encode the values in a standard eight bit format. Using
the normalization function the ratio value 1 is assigned the brightness value 128 and
the ratio values between 1/255 and 1 are assigned values between 1 and 128 by the
function:
BVi , j , n = Int.[(BVi , j ,r × 127) + 1].
Remote Sensing 69

Fig. 12 Band ratio image showing mineralization in lithology

Ratio values from 1 to 255 are assigned values within a range of 128–255 by the
function

BVi , j , n = Int. ⎡⎣(128 + BVi , j ,r 2)⎤⎦ .

Three pairs of ratio images can be co-registered (aligned) and projected as color
composites. In individual ratio images and in these composites, certain ground
features tend to be highlighted, based on unusual or anomalous ratio values. For
example, an ore deposit may be weathered or altered so that a diagnostic surface
staining, called gossan, develops. This stain consists of hydrated iron oxide (rust)
that is normally yellow-brown. In band 3, this material reflects strongly in the red
but it is apt to be dark in band 4. The ratio quotient values for this situation tend,
therefore, to exceed 2–3, giving rise to a bright spot pattern in a 3/4 image (Fig. 12).
The yellows and reds in this composite (Fig. 12) denote areas of rock alteration and
mineralization.

6 Classification

6.1 Image Classification

The computer aided classification of the remotely sensed data is based on the
concept that a pixel is characterized by its spectral signatures, which vary in diffe-
rent wave bands. The spectral signatures of themes are assigned a gray level or
color which relates the Dn values to thematic information. The information extraction
process which analyzes the spectral signatures and assigns the pixel to thematic
categories based on similarity of Dn values is referred to as classification. There are
two types of classification:
• Unsupervised classification
• Supervised classification
70 D. Al Ajmi and Saif ud din

6.2 Unsupervised Classification

In this classification numerical operations are performed that search for natural
groupings of the spectral properties of pixels as examined in an image. The com-
puter selects the mean class and covariance matrices to be used in classification.
Once the data is classified, the classified data are assigned to some natural and
spectral classes and the spectral classes are converted to information classes of
interest. Some of the clusters are meaningless as they represent mixed classes of
earth surface materials. The unsupervised classification attempts to cluster the Dn
values of the scene into natural boundaries using numerical operations.
Unsupervised classification operates on the color composite made from bands 2,
3, and 4 specifying just six clusters (Fig. 13).
The light buff colors associate with the marine waters but are also found in the
mountains where shadows are evident in the individual band and color composite
images. Red occurs where there is some heavy vegetation. Dark olive is found
almost exclusively in the ocean against the beach. The orange, green, and blue
colors have less discrete associations.
The image in Fig. 14 shows bands 2, 3, and 4, in which 15 clusters are set up; a
different color scheme is chosen (Fig. 14).
In this image many individual areas represented by clusters do not appear to
correlate well. Unfortunately, what is happening is a rather artificial subdivision of
spectral responses from small segments of the surface. Another composite, bands
4, 7, and 1, shows a new classification with the same problems as the first, although
sediment variation in the ocean is better discriminated (Fig. 15).

Fig. 13 Imagery depicting unsupervised classification

Remote Sensing 71

Cluster 1
Cluster 2
Cluster 3
Cluster 4
Cluster 5
Cluster 6
Cluster 7
Cluster 8
Cluster 9
Cluster 10
Cluster 11
Cluster 12
Cluster 13
Cluster 14
Cluster 15

Fig. 14 An unsupervised classification image with 15 classes

Cluster 1
Cluster 2
Cluster 3
Cluster 4
Cluster 5
Cluster 6
Cluster 7
Cluster 8
Cluster 9
Cluster 10
Cluster 11
Cluster 12
Cluster 13
Cluster 14
Cluster 15

Fig. 15 An unsupervised classification image

72 D. Al Ajmi and Saif ud din

The unsupervised classification is too much of a generalization and the clusters

only roughly match some of the actual classes. Its value is mainly as a guide to
the spectral content of a scene to aid in making a preliminary interpretation
prior to conducting the much more powerful supervised classification proce-
dures. There are several types of unsupervised classification which are
explained briefly below:

6.2.1 Sequential Clustering

This clustering algorithm operates in a two-pass mode. In the first pass, the pro-
gram reads through the data set and sequentially builds clusters (group of points in
spectral space). There is a mean vector associated with each cluster. In the second
pass, a minimum distance classification to mean vector algorithm is applied, pixel
wise, where each pixel is assigned to one mean vector created in pass 1.

6.2.2 Pass 1: Cluster Building

During the first pass, the analyst may be required to supply four types of
information:
1. R, a radius in spectral space used to determine when a new cluster should be
formed.
2. C, a spectral space distance parameter used when merging clusters.
3. N, the number of pixels to be evaluated between each merging of the clusters.
4. Cmax, the maximum number of clusters to be identified by the algorithm.
These can be set top default values, if no initial human interaction is desired.
In the sequential clustering the data set are evaluated sequentially from left to
right (line 1, column 1). The brightness value associated with each pixel in the
image represent the mean data vector of cluster, it is an n-dimensional mean vector
where n represents the number of bands. If the spectral distance (D) between two
clusters is greater than R, the mean data vector of cluster 1 becomes the average of
first and second pixel brightness values. The spectral distance (D) is computed
using Pythagorean theorem. If the distance between two clusters is less than D then
the two clusters are merged together.
Pass 2: Assignment of pixels to one of the Cmax Clusters using Minimum-Distance
Classification Logic
The final cluster mean data vectors are used in a minimum-distance to means
classification algorithm to classify all the pixels in the image to one of the Cmax
cluster. It is necessary to evaluate the location of the clusters in the image, label and
combine. It is usually necessary to combine some clusters.
Cluster labeling is performed by interactively displaying all the pixels assigned
to an individual cluster, making it possible to identify their location and spatial
association with other clusters. This interactive visual analysis in conjunction with
Remote Sensing 73

the information provided in the scatter plot, allows one to group clusters into infor-
mation classes.

6.2.3 Statistical Clustering

This method of unsupervised classification takes into account the homogeneity

of neighboring groups of pixels, instead of considering individual pixels equally.
The algorithm uses 3×3 sets of contiguous pixels that have similar measurement
vectors. The assumption behind this algorithm is that contiguous, homogenous
pixels usually indicate a spatial pattern within the data that is worth classifying. The
process consists of two parts:
• Homogeneity within the window of pixels being considered;
• Cluster merging and deletion.

6.2.4 Homogeneity Parameters

Windows of nine pixels (3×3) are tested for homogeneity, from the upper left corner
of the data, moving one window at a time, so that the windows do not overlap. A
skip factor may be specified so that every xth window across and every yth window
down is tested. The mean and standard deviation are calculated for the nine pixels
in each band. These values are compared to the values entered.
• L, a lower bound for the standard deviation. This value is usually small, but not
equal to zero. Its primary purpose is to prevent clusters with a standard deviation of
zero in any band. A cluster with a standard deviation of zero will also have a covari-
ance of zero with any band, causing zeros to appear in the covariance matrix.
• U, is an upper bound for the standard deviation. The higher the standard devia-
tion is for one band, the less homogenous the data is in that band. Therefore,
the upper bound is a ceiling for the amount of variation within a window.
• The coefficient of variation V is an alternative test for homogeneity, based on the
mean of the cluster.

6.2.5 Isodata Classification

Isodata stands for iterative self-organizing data analysis technique. It is iterative as

it repeatedly performs an entire classification and recalculates statistics. Self-
organizing refers to the way in which it locates a cluster with minimum user input.
The Isodata classification uses minimum spectral distance to assign a cluster for
each pixel, it will begin with a specified number of arbitrary cluster means and then
it processes repetitively, so that those arbitrary means will shift to the mean of the
cluster in the data. Clusters with large variances are split and clusters with smaller
variances are merged.
74 D. Al Ajmi and Saif ud din

6.3 Supervised Classification

In supervised classification the analyst has to specify all the parameters. The steps
that are to be taken while attempting supervised classification are:
• An appropriate classification scheme must be adopted.
• Representative training sites may be selected.
• Statistics must be extracted from the training site spectral data.
• The statistics are analyzed to select the appropriate features (band) to be used in
the classification process.
• Select the appropriate classification algorithm.
• Classify the imagery into m classes.
• Statistically evaluate the classification accuracy.

6.3.1 The Classification Scheme

The categories of interest must be carefully selected and defined to successfully

perform digital image classification. It is essential to realize the fundamental differ-
ence between the information classes (defined by the analyst) and spectral classes
(inherent of sensor). The major point of difference between various classification
schemes is their emphasis and ability to convert spectral classes into information
classes of remote sensing data.

6.3.2 Training

The aim of training is to obtain sets of spectral data that can be used to determine
decision rules for the classification of each pixel in the whole image data set. The
training data for each class must be representative of all data for that class. Each
training site consists of many pixels, conventionally it is taken that if there are n
number of bands the number of pixels in each band is n+1. The mean, standard
deviation, variance, minimum value, maximum value, variance-covariance matrix
and correlation matrix for training classes are calculated, which represent the fun-
damental information on the spectral characteristics of all classes. Since for selec-
tion of appropriate bands only this information is not enough, thus feature selection
is used. The training sites are presented on true color map bands 1, 2, and 3 for 13
classes (Fig. 16).

6.3.3 Feature Selection

Feature selection is the process of discriminating each class of interest and to deter-
mine the bands in which a particular class is highlighted. Feature selection involves
both statistical and/or graphical analysis to determine the degree of seperatability
Remote Sensing 75

seawater
sediment1
sediment2
Baysediment
marsh
wavesurf
sand
urban1
urban2
sunlits1
shadows1
scrublands
grass
fields
trees
cleared

Fig. 16 Selection of training data from satellite image

between classes. Combinations of bands are normally ranked according to their

ability to discriminate each class from all others using n bands at a time (Table1).
We can deduce from this table that most of the signatures have combinations of
Dn values that allow us to distinguish one from another, depending on the actual
standard deviations (not shown). Two classes, urban 1 and cleared (ground), are
quite similar in the first four bands but apparently are different enough in bands 5
and 7 to suppose that they are separable. The range of variations in the thermal band
6 is much smaller than in other bands, suggesting its limitation as an efficient
separator.
On the basis of the above training sets are prepared and the image can be clas-
sified; there are several types of classification methods.

6.4 Parallelepiped Classification

The classification is based on simple Boolean logic. Training data in n spectral

bands are used in performing the classification. Brightness values from each pixel
of the multi-spectral imagery are used to produce an n-dimensional mean vector,
76 D. Al Ajmi and Saif ud din

Table 1 Table of band means and sample size for each class training set
Class/BAND 1 2 3 4 5 6 (TH) 7 No. of pixels
1. Seawater 57.4 16.0 12.0 5.6 3.4 112.0 1.5 2,433
2. Sediments1 62.2 19.6 13.5 5.6 3.5 112.2 1.6 681
3. Sediments2 69.8 25.3 18.8 6.3 3.5 112.2 1.5 405
4. Bay sediment 59.6 20.2 16.9 6.0 3.4 111.9 1.6 598
5. Marsh 61.6 22.8 27.2 42.0 37.3 117.9 14.9 861
6. Waves surf 189.5 88.0 100.9 56.3 22.3 111.9 6.4 1,001
7. Sand 90.6 41.8 54.2 43.9 86.3 121.3 52.8 812
8. Urban1 77.9 32.3 39.3 37.5 53.9 123.5 29.6 747
9. Urban2 68.0 27.0 32.7 36.3 52.9 125.7 27.7 2,256
10. Sun slope 75.9 31.7 40.8 43.5 107.2 126.5 51.4 5,476
11. Shade slope 51.8 15.6 13.8 15.6 14.0 109.8 5.6 976
12. Scrublands 66.0 24.8 29.0 27.5 58.4 114.3 29.4 1,085
13. Grass 67.9 27.6 32.0 49.9 89.2 117.4 39.3 590
14. Fields 59.9 22.7 22.6 54.5 46.6 115.8 18.3 259
15. Trees 55.8 19.6 20.2 35.7 42.0 108.8 16.6 2,048
16. Cleared 73.7 30.5 39.2 37.1 88.4 127.9 45.2 309

Mc. Where Mc=(mc1, mc2, mc3,…,mcm) with mck, being the mean value of the training
data obtained for class c in band k out of m possible classes. Sck is the standard
deviation of the training class c in band k out of m possible classes. Using one
standard deviation threshold, a parallelepiped algorithm decides BVijk is in class
c, if

m ck − Sck ≤ BVijk ≤ m ck + Sck ,

where c=1,2,3,…,m number of classes, k=1,2,3,…,n number of bands; thus the low
and high decision boundaries are defined as

Lowck = m ck − Sck .

Highck=mck+Sck the parallelepiped algorithm becomes

Lowck ≤ BVijk ≤ High ck .

These decision boundaries form an n-dimensional parallelepiped in future space. If

the pixel value lies above the lower threshold and below the high threshold for all n
bands evaluated it is assigned to that class. When an unknown pixel does not satisfy
any of the Boolean logic criteria, it is assigned to an unclassified category. Increasing
the size of thresholds to classify the unclassified category, would increase the size of
the parallelepipeds and introduce a significant amount of overlap among many paral-
lelepipeds resulting in classification error. Perhaps these pixels belong to the class that
Remote Sensing 77

has not been trained. The parallelepiped algorithm is a computationally efficient

method, unfortunately because some parallelepipeds overlap, it is possible that an
unknown candidate pixel might satisfy the criteria of more than one class, in such
cases the pixel is assigned to the first class for which it meets all criteria [3]. The
parallelepiped classification in multi-band data sometimes resulted in loss of thematic
classes due to similar spectral expression of more than one class.

6.5 The Minimum Distance to Mean Classification

This is a simple and commonly used classification algorithm and the classification
accuracies are comparable to any other classification. The user in this classification
has to provide the mean vectors for each class in each band mck, from the training
sets. In this classification distance of each mean vector mck is calculated for each
unknown pixel (BVijk), this distance is calculated using Euclidian distance based on
Pythagorean theorem or “around the block” distance measures. The computation of
the Euclidian distance from point a (x,y) to the mean of class 1 (mx, my) measured in
band k and l, relies on equation:

Distance = √ (BVijk − m ck )2 + (BVijl − m cl )2 ,

where mck and mcl represent the mean vectors of class c measured in bands k and l.
The subscript for class c is incremented from 1 to n, by calculating the Euclidian
distance from point a to the mean of all the classes, it is possible to determine which
distance is shortest. In minimum distance a threshold can be assigned from the class
means, beyond which a pixel will not be assigned to a category even though, it is
nearest to the mean of that category. When more than two bands are evaluated in a
classification, it is possible to extend the logic of computing the distance between
just two points in n-space using equation.
n
Distance = √ ∑ (a1 − b1 )2 .
i =1

Each unknown pixel is then placed in the class closest to the mean vector in this
band space. For this classified image there were 16 gray levels, each representing a
class, to which a color is assigned. This minimum distance classification has all
seven TM bands, including the thermal (Fig. 17).

6.6 The Maximum Likelihood Classification

This classification assigns each pixel having feature x to the class c whose units
are most probable to have given rise to feature vector x. It assumes that training
data statistics for each class in each band are Gaussian in nature (Uni-modal).
78 D. Al Ajmi and Saif ud din

seawater
sediment1
sediment2
Baysediment
marsh
wavesurf
sand
urban1
urban2
sunlits1
shadows1
scrublands
grass
fields
trees
cleared

Fig. 17 Supervised classification of satellite image

Bi-modal or tri-modal histograms in a single band are not ideal for max-like
classification, in such cases, individual modes probably represent individual
classes that should be trained upon individually and labeled as separate classes,
thus producing unimodal data.
Maximum likelihood classification makes use of the mean measurement vector,
Mc, for each class and covariance matrix for class c for bands k through i, Vc.
pc≥pi, where i=1,2,3,…m possible classes

pc = {−0.5 log[det (Vc )]} − {0.5 ( X − M c )T (Vc−1 ) ( X − M c )},

where det (Vc) is the determinant of the covariance matrix (Vc).

To classify the measurement vector x of an unknown pixel into a class, the maxi-
mum likelihood decision rule computes the value pc for each class. Then it assigns
the pixel to the class that has the maximum value.
In this image 16 classes are identified. These classes are identical to the previous
ones recorded in the minimum distance image. In both instances, the Sediment class
has been subdivided into three levels and two Urban classes are attempted, to
account for visual differences between them (Fig. 18).
Remote Sensing 79

seawater
sediment1
sediment2
baysediment
marsh
wavesurf
sand
urban1
urban2
sunlits1
shadows1
scrublands
grass
fields
trees
cleared

Fig. 18 Supervised classification images with additional classes in sediments and urban area

7 Monitoring Applications of Environmental

Consequences of War

7.1 Application I – Hydrocarbon Pollution consequent

to 1991 Gulf War

A unique environmental catastrophe affected Kuwait during the 1991 Gulf War,
which severely impacted environmental conditions. The most conspicuous among
these changes were development of oil lakes and oil polluted surfaces. The oil lakes
were formed due to gushing of oil from the free flowing oil wells. Thick hydrocar-
bon deposits covered large areas in the form of tarmats and tarcrete over the desert
surface from the burning oil. Studies post the 1991 Gulf War estimated around 300
oil lakes covering an area of 36km2. These interpretations were largely based on
remote sensing studies with selective ground checks. The temporal variation in oil
lakes has been studied in detail by Kwarteng and Al-Amji [19]; their study shows
a decrease in spatial coverage in recent years.
Monitoring near surface existence of hydrocarbon polluted surfaces in the
Burgan Oil field area was attempted by Saif ud din et al. [20, 21]. The researchers
used the land surface temperature (LST) as an indicator for hydrocarbon pollution.
80 D. Al Ajmi and Saif ud din

The LST and total petroleum hydrocarbon (TPH) are positively correlated. The
methodology followed is based on the fact that most of the satellites carry a ther-
mal infra-red band which can be used for LST estimation. The spatio-temporal
variation in the thermodynamic properties of surface material has been mapped in
order to identify hydrocarbon polluted surfaces using Landsat TM data. Emissivity
is a strong indicator of compositional variation in silicate minerals which make up
the bulk of the earth’s surface material. Emissivity affects the apparent tempera-
ture due to changes in the thermal properties of materials (conductivity, density,
capacity, and inertia). There are several algorithms proposed to estimate LST from
remotely sensed data. The most common of these are mono-window and split
window methods [22–25], the latter was used initially to estimate sea surface
temperature.
The land surface temperature measurement is a complicated task due to high
spatio-temporal variation of surface emissivities and atmospheric water vapor as
both affect the thermal radiance reaching the sensor [26]. Landsat TM band 6 data
with a wavelength range of 10.45–12.50 mm can be utilized to estimate land surface
temperature. The Landsat TM has a single thermal band, therefore the split window
method [27], and temperature–emissivity method [28], could not be applied. For
measurement of LST from Landsat TM data, the mono-window algorithm, single
channel algorithm, and radiative transfer equations, can be used [23, 25, 29, 30].
The algorithm proposed by Richter [23] is simple and accurate. The logic for LST
retrieval from satellite data is based on the fact that ground emissivity is known.
LST can therefore be calculated by accounting for atmospheric correction.
Computation of atmospheric correction is a complex process. To simplify the proc-
ess, however, upward atmospheric thermal radiance and the reflected atmospheric
radiance are subtracted from the observed radiance at satellite to compute bright-
ness temperature at ground level. The standard atmospheric profiles provided in the
ATCOR program of PCI, which is based on LOWTRAN can be used to get the
atmospheric radiance for surface temperature estimation using band 6 data from
Landsat TM.
The mono-window algorithm proposed by Qin et al. [25] is based on the thermal
radiance transfer equation to calculate LST. It utilizes transmittance and mean
atmospheric temperature to estimate LST. The LST estimation is done considering
the fact that brightness temperature at satellite can be computed by estimation of
radiance from Dn value and conversion of radiance into brightness temperature.
The radiance calculation from Dn of TM data utilizes an equation developed by
Markham and Barker [31], shown below.

(
L( l ) = Lmin( l ) + Lmax( l ) + Lmin( l ) )Q
DN
Qmax
,

where L(l), spectral radiance received at sensor; Qmax, maximum Dn value, i.e., 255;
QDn, Dn value of the selected pixel, i.e., it can be anywhere between 0 and 255;
Lmin(l) and Lmax(l) are minimum and maximum spectral radiance for QDn=0 and
Qmax=255.
Remote Sensing 81

The ATCOR Program of PCI uses the constants defined by Schneider and Mauser
[32], where the average wavelength for Landsat TM band 6 is taken as 11.475 μm,
Lmin(l) = 0.1238 mWcm−2 sr−1 μm−1 for Dn value 0 and Lmax(l)=1.56mWcm−2sr−1 μm−1
for Dn value 255. Substituting these values in the equation above

(
L( l ) = 0.1238 1.56 − 0.1238
255 )Q
dn

L( l ) = 0.1238 + 0.005632156Qdn .

The brightness temperature correction for true LST is based on the radiative trans-
fer equation [33]. The concept accounts for thermal emittance from an object in
accordance with the blackbody theory, which states:
C1
B( l ) (T ) =
( )
C2
,
5
l e lT
−1

where B(l)(T) is spectral radiance of a blackbody measured in Wm−2sr−1 μm−1; l is

wavelength in meters; C1 and C2 are spectral constants; C1=1.19104356×10−16 Wm2;
C2=1.4387685×104 μmoK; T=temperature oK.
This concept of a blackbody is theoretical. In real-life situations, materials eve-
rywhere do not behave as blackbodies because of their emissivity.
Emissivity = Fr / Fb,

where Fr, radiant flux exiting a real-world body; Fb, radiant flux exiting a black
body.
Since the materials in the study area do not behave as blackbodies, emissivities
are to be considered which will account for emittance by atmospheric absorption
and reflection.
The LST estimates take into account the upward and downward atmospheric
radiances. The former is greater than the latter. Their difference in clear sky is
within 5°C [25]. In the present study, the emissivity for the sand in a desert environ-
ment is assumed to be 0.76, for carbonate lithologies in the study area it is 0.92, and
for heavy oil it is 0.97. The approximation is taken as:

B6 (Ta∞ ) = B6 (Ta ).

To get LST from Ts Planck’s radiance function is linearized [25].

L6 = a + bT6 ,

where a and b are coefficients valid for a temperature range of 0–70°C.

a = −67.355351
b = 0.45860
The ATCOR program used for estimation of the LST assumes two surface tempera-
tures Ts1 and Ts2 with their corresponding black body temperatures Tb1 and Tb2.
82 D. Al Ajmi and Saif ud din

The atmospheric correction and elevation functions are included in ATCOR, which
provide a very good approximation. The temperature images are derived from
ATCORT1.
The digital image processing of Landsat TM band 6 for land surface temperature
measurement carried out for post-1991 Gulf War images show a remarkable pattern
coincident with the surface hydrocarbon disposition in the area expressed as vari-
able LST within the image. The area shows an increase in the LST as compared to
a pre-war image of 1989.
Intra-image LST variations depict oil polluted surfaces since the emissivity
and thermodynamic behavior of oil is different from that of soil or water. The
1989 image shows LST differences due to different emissive properties of the
lithology and sand. The area exhibits LST variations in the sand cover in the post-
1991 Gulf War images, due to hydrocarbon polluted surfaces formed during the
1991 Gulf War. The persistent spatial intra-image LST pattern correlates with the
oil lakes, tarmats, and soot fallout. Selective ground truth has been carried out by
Kuwait Institute for Scientific Research [34–38] (Figs. 19–28).
The intra-image LST variations were verified from the collated ground truth data
to verify the temperature anomalies in the images (Fig. 29).
The intra-image LST variations correlate very well with the hydrocarbon pol-
luted sites in the study area. The algorithm applied to study intra-image spatial vari-
ations in LST shows temperature variations in 1991, 1995, 1998, and 2000 in the
Burgan Oil Field of Kuwait. These figures show a similar LST pattern during
successive years even though they were obtained at different times of the year. The
LST pattern correlates with the hydrocarbon polluted surfaces in the study area. A
sketch map is included to show different pollution zones and location of some field
photographs is also marked on the LST map. The temporal changes have resulted
in reduction of spatial extent and pollution intensity of these zones. The reduction
in intensity is primarily due to weathering of the hydrocarbon, leaching during rains
and downward migration in summers when surface temperature is extremely high
(exceeds 50°C).
Intra-scene LST values are directly proportional to pollution intensity. The
higher temperatures in the LST maps indicate relatively higher pollution. It is
clear from Figs.15–18 that the intensity of pollution and spatial extent has
reduced in successive years. This decrease is probably due to several reasons:
increase in depth of pollutants; bioremediation; sludge recovery; leaching; and
chemical alterations.
Some of the “oil lakes” and tarcrete deposits are covered by a thin veneer of
sand which has obscured them. The depth of burial has increased over the
years; however, the thickness of sand sheet has not increased in recent years in
the study area [39]. High temperatures in summer months, imparts mobility to
the tar and soot deposits at the surface, and winter rains leach out hydrocarbon
fractions which percolate down with the water [40]. Consequently, the depth of
pollution observed at several places, varies from a surface veneer to a maximum
depth of 4μ.
Remote Sensing 83

3230000
W E

S
3225000

Ahmadi

3220000

Quarry
3215000

3210000

3205000

3200000

3195000

3190000

775000 780000 785000 790000 795000 800000

0 5 Kilometer

Fig. 19 Landsat image over Burgan Oil field (1989)

The analyses of pre- and post-1991 Gulf War images of the study area show a
very elegant co-relation of LST with the spatial locations of oil lakes, tarcrete, tar-
mat, and soot in the Burgan oil field. Higher intra-image LST corresponds to higher
84 D. Al Ajmi and Saif ud din

Legend
N Temperature in
Celsius
14 - 16
3230000
W E
16 - 18
18 - 20
S
3225000 20 - 22
22 - 24
Ahmadi
24 - 26
3220000 26 - 28

28 - 30
Quarry 30 - 32
3215000
32 - 34
34 - 36

3210000
36 - 38
38 - 40

40 - 42
3205000 42 - 44

3200000

0 5 Kilometer
3195000

3190000

775000 780000 785000 790000 795000 800000

Fig. 20 LST image over Burgan Oil field (1989)

total petroleum hydrocarbon (TPH) and/or near surface burial. Very high values are
also observed over the flaring chimneys of the oil separators, while lower LST
correlates well with lower TPH and/or deep burial and water pits. The TPH concen-
tration in the uncontaminated soil varies from 7 to 20 μgg−1 whereas the “oil lakes”
shows a TPH concentration of around 25,000 μgg−1. These “oil lakes” correspond
Remote Sensing 85

W E
3230000
S

3225000

Ahmadi
3220000

3215000

Oil Lakes

3210000 1
2

3205000

3200000

3195000

3190000

775000 780000 785000 790000 795000 800000

0 5 Kilometer

Fig. 21 Landsat image over Burgan Oil field (1991)

to the high temperature zones in LST maps. The tarmats are solidified thick sludge
deposits with TPH concentration of 12,000–19,000 μgg−1. They correspond to a
slightly lower temperature than “oil lakes.” There exists a good correlation between
the hydrocarbon concentrations and LST in all images.
86 D. Al Ajmi and Saif ud din

N Legend
Temperature in
W E
Celsius
3230000 14 - 16
S 16 - 18
18 - 20
3225000 20 - 22
22 - 24
Ahmadi
24 - 26
3220000
26 - 28
28 - 30

3215000 30 - 32

Oil Lakes 32 - 34
34 - 36
3210000 1
2 36 - 38
38 - 40

3205000 40 - 42
42 - 44

3200000

3
3195000 0 5 Kilometer

3190000

775000 780000 785000 790000 795000 800000

Fig. 22 LST image over Burgan Oil field (1989)

The LST accuracy is ± 0.9°C to 1.1°C. However, in this study the accuracy of
the actual LST estimation is not a primary concern since the relative LST patterns
are sufficient to identify different hydrocarbon polluted surfaces. The actual LST
measurement and calibration is part of an ongoing research.
The methodology developed can be adopted for detecting oil polluted surfaces
on land and as spills in oceans and near the ports [20]. The technology can be used
Remote Sensing 87

3230000 W E

S
3225000
Ahmadi

Water pit
3220000

3215000

3210000

Oil Lakes
3205000

3200000

3195000

3190000
Road

775000 780000 785000 790000 795000 800000

0 5 Kilometer

Fig. 23 Landsat image over Burgan Oil field (1995)

88 D. Al Ajmi and Saif ud din

N Legend
Temperature in
W E Celsius
3230000 14 - 16
S
16 - 18
18 - 20
3225000 20 - 22
Ahmadi 22 - 24
24 - 26
Water pits
3220000 26 - 28
28 - 30
30 - 32
3215000
32 - 34
34 - 36
3210000 36 - 38
38 - 40
Oil lakes 40 - 42
3205000
42 - 44

3200000

0 5 Kilometer
3195000

3190000
Road

775000 780000 785000 790000 795000 800000

Fig. 24 LST image over Burgan Oil field (1995)

for quick monitoring of spatially large areas, using similar algorithms for other
satellite data, where the temporal coverage is quick. An important application can
be in near real time monitoring of oil spills. The effect of darkness will not affect
the LST observations so waste dumping during night-time can be easily identified
for pollution control and environmental management (Fig.30).
Remote Sensing 89

3230000 W E

S
3225000
Ahmadi

3220000

Quarry

3215000

3210000

3205000

3200000 Oil lakes

3195000

3190000

775000 780000 785000 790000 795000 800000

0 5 Kilometer

Fig. 25 Landsat image over Burgan Oil field (1998)

90 D. Al Ajmi and Saif ud din

Legend
N Temperature in
Celsius
3230000 W E 14 - 16
16 - 18
S 18 - 20
3225000 Ahmadi
20 - 22
22 - 24
24 - 26
3220000
Quarry 26 - 28
28 - 30
3215000
30 - 32

32 - 34

3210000 34 - 36
36 - 38
38 - 40
3205000
40 - 42
42 - 44

3200000
Oil lakes

3195000
0 5 Kilometer

3190000

775000 780000 785000 790000 795000 800000

Fig. 26 LST image over Burgan Oil field (1998)

7.2 Application II – Precipitation Estimate Using Tropical

Rainfall Measuring Mission

In the post-war scenario the precipitation patterns changed. The total rainfall meas-
uring mission data has been calibrated over Kuwait to have continuous spatial and
temporal precipitation coverage over Kuwait. Precipitation retrievals from remote
Remote Sensing 91

3230000 W E

S
3225000
Water pit Ahmadi

3220000

Quarry

3215000

3210000 Oil lakes

3205000

3200000

3195000

3190000

775000 780000 785000 790000 795000 800000

0 5 Kilometer

Fig. 27 Landsat image over Burgan Oil field (2000)

92 D. Al Ajmi and Saif ud din

N Legend
Temperature in
Celsius
3230000 W E 14 - 16
16 - 18
S 18 - 20
3225000 20 - 22
Ahmadi
Water pits 22 - 24
24 - 26
3220000
26 - 28
Quarry
28 - 30
3215000 30 - 32

32 - 34
Oil lakes 34 - 36
3210000
36 - 38
38 - 40
3205000 40 - 42
42 - 44

3200000

3195000
0 5 Kilometer

3190000

775000 780000 785000 790000 795000 800000

Fig. 28 LST image over Burgan Oil field (2000)

sensing satellite data are an acceptable alternative to the use of rain gauges for
global and regional climatic/hydrological models. The tropical rainfall measuring
mission (TRMM) is a low altitude satellite, equipped with an active precipitation
radar (PR) along with a multi-channel passive TRMM microwave imager (TMI),
visible and infrared scanners (VRIS), earth radiant energy flux system, and light-
ning imaging sensor. TRMM provides better spatio-temporal precipitation cover-
age on global and regional scales, as it provides a wide variety of data.
Remote Sensing 93

3230000 W E

S
3225000

3220000

3215000

3210000

3205000

3200000

3195000

3190000

775000 780000 785000 790000 795000 800000

0 5 Kilometer

Fig. 29 Spatial distribution of soot, tarmats, tarcrete and oil lakes in southern Kuwait

The TRMM precipitation retrievals are based on the emission of microwave

radiation from raindrops. The PR is the first space-based radar to measure precipita-
tion [41–43]. It operates at 13.8 GHz frequency with a 2.17 cm wavelength. The
instrument is capable of detecting echoes of >17 dbz. However, due to the strong
reflection of the earth’s surface, PR measures precipitation rates close to the ground
surface, thus affecting the reflectivity profile above the surface [43].
In developing countries, particularly in arid countries, there is a paucity of con-
tinuous rain gauge precipitation measurements from which time series analyses can
be carried out if required. An attempt is being made to redress this balance, in
Kuwait. A technique to resample the monthly TRMM precipitation rate on 0.25° is
described for calculating precipitation at any spatial location with limited or no
ground control [44]. The technique has been calibrated by testing modeled data
94 D. Al Ajmi and Saif ud din

3230000 W E Legend
Concentration in microgram/gram
S (Temperature in degree Celcius)
3225000
20 - 50 (19 - 21)
136 - 1004 (21 - 23)
3220000
4114 - 6827 (25 - 27)
8650 - 9440 (27 - 29)
3215000
12400 - 18380 (31 - 33)

25390 - 26300 (35 - 37)

3210000

3205000 0 5 Kilometer

3200000

3195000

3190000

775000 780000 785000 790000 795000 800000

Fig. 30 Correlation between surface temperature and hydrocarbon pollution

with actual rain gauge measurements at eight different locations in Kuwait for
which sufficient ground measurements are available over a 96-month period, for the
period between 1st January, 1998 and 31st December, 2005.
The Kuwait Institute for Scientific Research operates eight weather stations for
gathering baseline metrological data over Kuwait. The stations are equipped with
MetOne 370 rain gauges that are operated on a tipping bucket principle. This allows
accurate repeatable measurement with minimum maintenance. The data is recorded
every 10 minutes and is stored in attached data loggers. TRMM dataset 3B43 V6 is
used to obtain satellite precipitation estimates. The 3B43 is executed once every
calendar month and gives a single best estimate for precipitation rates and the root
mean square precipitation error estimate, by combining three-hourly integrated
high quality data, IR estimate (3B42) with the monthly accumulated climate assess-
ment monitoring system (CAMS) or global precipitation climatology center
(GPCC) rain gauge analyses (3S45). This dataset is corrected on both a global and
regional scale.
Remote Sensing 95

In this application, an attempt has been made to validate weighted bi-linear

interpolation for precipitation retrievals over a point of interest. The source dataset
is the TRMM 3B43 V6 data with a 0.25 degree grid. This grid is a bit coarse for
local interpretation. Instead of using the grid value for a particular location, or the
averaging of adjacent pixel values on a 0.25 degree grid, we found that the bi-linear
weighted interpolation which is presented here will be more realistic. The TRMM
3B43 dataset was used for bilinear weighted interpolation to resample the data for
a specific location [45], which in the present study is the spatial location of the rain
gauge. The logic that has been adopted is to average four adjacent pixels for a location
of interest [46]. The location (X,Y) has been kept in a 2×2 grid in such a manner
that it occupies a central position in the grid. An illustration is included (Fig. 31).
In the bi-linear interpolation, a simplistic rationale is followed, i.e., for calculat-
ing the pixel value of a particular position (X,Y), four adjacent pixel values are used
[21]. The closer the pixel is to the position (X,Y), the more influence (weight) it will
carry. The method is not merely a falling function of distance from the pixel. Rather
it considers a weighted approach based on its spatial locations in a two-dimensional
space [47]. The derivation of bi-linear interpolation weights can be expressed as
follows [48, 49]:
X = S x X1 + (1 − S x ) X 2 ,
X − X1
Sx = ,
X 2 − X1
Y = S yY1 + (1 − S y )Y2 ,
Y − Y1
Sy = ,
Y2 − Y1

X1Y1 X2Y1
X

(X,Y) Y

X1Y2 X2Y2

Fig. 31 Grid considered for estimation of point values

96 D. Al Ajmi and Saif ud din

where 0 ≤ S ≤ 1 ⇔ X1 ≤ X ≤ X 2 .

Considering the above, the actual weight at any point (X,Y) in a two-dimensional
space can be computed [49] as:

I ( X , Y ) = (1 − S x )(1 − S y )I ( X1Y1 ) + S x (1 − S y )I ( X 2Y1 ) + (1 − S x )S y I ( X1Y2 ) + S x S y I ( X 2Y2 ),

where I=actual pixel value.

The comparison of the satellite precipitation following the bilinear weighted
interpolation and the gauge measurement are presented below at eight different
spatial locations. The measured data for the rain gauges was not continuous, and at
certain times was not recorded due to instrument malfunction and/or calibration
errors. However, the TRMM data is continuous (Figs. 32–39).
The bi-linear weighted interpolation has been applied to the entire time series of
TRMM dataset for 96 months to get calculated data over each station even when
the measurement was not made by gauge. It is obvious, from the plots, that the
measured and computed precipitation data are very highly correlated. This suggests
that the methodology can be a reliable alternative to assessing precipitation in areas
where rain gauge measurements are unavailable. Normally, the 3B43 data set is
corrected using ground truth measurements, but their validation sites are spaced far
apart with some even located on ocean buoys, which may not even account for the
precipitation variation on a local scale due to Urban Heat Island, variable emissivity
of surface material and surface topography.
Rain gauge measurements at all eight locations were compared to: (1) pixel values
for the grid within which rain gauges were cited, (2) an average 2×2 grid value around
a particular rain gauge location and (3) a bi-linearly weighted interpolation estimate

Taweel
100
90
80
Rainfall in Millimeter

70
60
50
40
30
20
10
0
1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93
Months
Rain Gauge TRMM Estimate

Fig. 32 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at Taweel
Remote Sensing 97

Salmi
120

100
Rainfall in Millimeter

80

60

40

20

0
1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93
Months
Rain Gauge TRMM Estimate

Fig. 33 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at Salmi

KISR
160

140

120
Rainfall in Millimeter

100

80

60

40

20

0
1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93
Months
Rain Gauge TRMM Estimate

Fig. 34 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at KISR

around a rain gauge of interest. It was found that the rain gauge measurements at all
eight locations compared extremely well with method (3) above. This suggests that
more accurate estimates of precipitation at any location are obtained by the use of the
bilinear weighted interpolation methodology. This study therefore validates the use of
98 D. Al Ajmi and Saif ud din

Azor
120

100
Rainfall in Millimeter

80

60

40

20

0
1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93
Months
Rain Gauge TRMM Estimate

Fig. 35 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at Azor

Mutla
120

100
Rainfall in Millimeter

80

60

40

20

0
1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93
Months
Rain Gauge TRMM Estimate

Fig. 36 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at Mutla
Remote Sensing 99

Subbiya
120

100
Rainfall in Millimeter

80

60

40

20

0
1 6 11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91 96
Months
Rain Gauge TRMM Estimate

Fig. 37 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at Subbiya

Haiman
90
80
Precipitation in Millimeter

70
60
50
40
30
20
10
0
1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93
Months
Rain Gauge TRMM Estimate

Fig. 38 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at Haiman

a bi-linear weighted interpolation method on a TRMM 3B43 V6 dataset. The com-

puted precipitation rates are slightly higher than measured values albeit within +2%
of the measured values. This is extremely good considering that rain gauges are
believed to underestimate precipitation of the order of 10–20% due to the wind speed,
precipitation rates, and sometimes the apparatus itself [50–54].
100 D. Al Ajmi and Saif ud din

Wafra
120

100
Rainfall in Millimeter

80

60

40

20

0
1 5 9 13 17 21 25 29 33 37 41 45 49 53 57 61 65 69 73 77 81 85 89 93
Months
Rain Gauge TRMM Estimate

Fig. 39 Correlation between precipitation estimates obtained from Tropical rainfall measuring
mission satellite and rain gauge at Wafra

Fig. 40 An oil logged farm

7.3 Application III – Feature Extraction Technique for Palm

Tree Census

Vegetation monitoring post the 1991 Gulf War is of prime importance. Lots of date
palm trees and farms have been extensively damaged (Fig. 40).
Remote sensing technology is being used to automatically classify the vegeta-
tion type for agriculture census. The Kuwait Institute for Scientific Institute took a
step in this direction by developing a methodology for automatic mapping of urban
treed areas. During the current phase it is mode directed towards mapping of date
Remote Sensing 101

palm trees. Palm trees are very common in Middle Eastern countries. They are of
significant environmental and commercial importance [55]. In recent decades the
Middle Eastern region has witnessed an extensive planting of date palm trees, both
in urban and agricultural areas. Millions of trees are estimated to have been planted
in these arid deserts. Extensive plantation in urban areas rarely gives a clue that
these are arid and hyper arid countries. Among these species most common are date
palm trees, which are seen planted along the roads, in front of houses, in parks, and
organized plantation in agricultural areas. However, there is limited knowledge of
actual tree counts and their exact spatial locations, which is a requirement for any
agricultural census.
Remote sensing data has been used for the identification of urban treed areas,
but with limited classification accuracies. These lower classification accuracies are
attributed to a variety of spectral and textural properties [56]. Medium resolution
satellites including LANDSAT, SPOT, and ASTER have been used in urban treed
classification but their spatial resolution permits only larger patches of treed areas
to be classified [56–58]. With the advancement in Satellite technology and availa-
bility of high spatial resolution images, it is now feasible to achieve higher classi-
fication accuracies in urban and agricultural areas. The present study is an attempt
to map the date palm trees in urban and agricultural areas of Kuwait. An accuracy
assessment is also made to compare results from maximum likelihood and the
Laplacian blob classifications of date palm trees within the test areas.
Similar studies for classifying and quantifying olive trees were taken up by
European Union Countries. The European Economic Committee (EEC) realized
the need to quantify the Olive plantation in 1997 and launched the OLISTAT
project in September 1997 to estimate the number of olive trees in France, Italy,
Spain, Portugal, and Greece [59]. The counting of trees is a classic example of
remote sensing applications in forestry. However, crown counting of trees is not an
easy or straightforward task as there are limitations of satellite data resolution as
well as problems related to the subjective nature of interpretation. Howard [60]
indicated that the capacity to distinguish different objects is governed by size of the
object relative to pixel.
The multi-spectral classification methods have provided reasonably good
results, but there is still room for further improvement in classification accuracy if
textural parameters are taken into account. It was believed that with the availability
of higher resolution satellite data, classification accuracies will improve, but
~100% accuracies are yet to be seen. Since the higher spectral resolutions increased
the intra-class separatability in an image, Marceau et al. [61] suggested that optimal
spatial resolution for the classification of temperate forest should be 10m as resolu-
tions finer than 10m increased the intra-class seperatability and decreased the
overall classification accuracy.
The inclusion of textural parameters in urban and forest area image classification
improved accuracies of image classification on high resolution images [56, 62].
Previous investigations have attempted the extraction of tree textures from high-
resolution satellite data using neural network, co-occurrence matrix, semi-variogram,
threshold based spatial clustering, local variance and local maximum filtering
102 D. Al Ajmi and Saif ud din

[63–68]. Some of these approaches worked well but the success rate was limited in
urban areas. The importance of urban area classification in this region is realized
since there is extensive date palm plantation in urban areas throughout Middle
Eastern countries. As these trees are planted along the roads, inside and outside
private properties and on road side pavements, the gathering of information of these
trees for agricultural census purposes is a cumbersome task to achieve, yet there is
no systematic database, to document their spatial information.
Researchers started to develop algorithms and application software to count
trees [69–71]. In this communication, selective filtering and Laplacian blob are
used for classifying date palm trees in both the urban and agricultural areas.
Two Quickbird scenes of April 2005 were selected over the study area. One of
the scenes is over an urban area in central Kuwait, while the other is from an agri-
cultural area in northern Kuwait. The quickbird data set used is panfused with a
spatial resolution of 0.6m (Fig. 41).
There are two steps involved in palm tree counting. The first step is the selective
smoothing by non-linear diffusion, which results in a sharp contrast among a
number of different features followed by Laplacian filtering.

Iraq Iran

Kuwait

46.6 46.8 47 47.2 47.4 47.6 47.8 48 48.2 48.4 48.6

Qatar 30.2 30.2
UAE
Saudi Arabia N

30 30
an
m

W E
O

29.8 S 29.8
n
me
Ye

29.6 29.6

29.4 29.4
Kuwait

29.2 29.2

29 29

28.8 Study Area 28.8

28.6 28.6
0 30 km
28.4 28.4
46.5 47 47.5 48 48.5

Fig. 41 Location of the two selected scenes

Remote Sensing 103

The second step involves the application of a non-linear parabolic equation [72]
to the two selected quickbird scenes. This equation allows selective enhancement
and smoothing in addition to simultaneously preventing the blurring of the edges.
The processing is quite effective for image classification in urban areas. The equa-
tion is stated as:
⎛ ∇I ⎞
∂ I ( x, y, t ) / ∂ t = g (Gx∇I ) ∇I div ⎜ ⎟,
⎝ ∇I ⎠
⎛ ∇I ⎞
where ∇I div ⎜ ⎟ diffuses the image I ( x, y) in the direction orthogonal to its
⎝ ∇I ⎠
gradient ∇I and not in all directions.

[g (Gx∇I ) is used for edge enhancements.

This anisotropic filtering is basically the statistical interpretation of the aniso-
tropic diffusion, but it takes slightly more time in processing and implementation
(Black and Sapiro 1998).
The second step is Laplacian filtering which can be suggested as an irreducible
differential invariant. It is expressed mathematically by equation
[∇2 I = Ixx + Iyy,

when a gray scale image defines a first-order derivative as equation ∇I = (Ix, Iy )

⎛ Ixx Ixy⎞
and a second-order derivative by Hessian matrix H I = ⎜ .
⎝ Iyx Iyy⎟⎠

When the second-order derivative is greater than zero and I x + I y = 0 , then the
2 2

point is referred to as an elliptic point due to its appearance. The sign of Laplacian
equation indicates the maxima and minima. The dark blob indicates the minima (A)
and the bright blob indicates the maxima (B) (Fig. 42).

I xx + I yy > 0.
I xx + I yy < 0.

The proposed methodology is morphometric and thus spatial resolution in palm

tree mapping is critically important. The use of 0.6-m spatial resolution data to map
palm trees with 3–4 m crown sizes and 3–8 m inter tree spacing was sufficient.
Similar methodology has been successfully used in the European Union countries
for the mapping of Olive trees. Application of this methodology to classify date
palm trees in arid urban and agriculture areas was employed using two Quickbird
scenes. The original quickbird scene with different species of trees and bushes with
similar spectral characteristics, including streets and cars and the other scene with
green houses, buildings and trees are used for classifying using maximum likeli-
hood classification. The results of this classification are not perfect. A simple
104 D. Al Ajmi and Saif ud din

1
0.5
−1 0
−0.5 −0.5
0
0.5 −1
1
Ixx + Iyy > 0
b

1
0.5
−1 0
−0.5 −0.5
0
0.5 −1
1
Ixx + Iyy < 0

Fig. 42 Laplacian blob (a) minima and (b) maxima

accuracy assessment was carried out using the random point method in an area of
1,000×1,000 pixels for images.
A number of 100 randomly selected points were assessed in either image to
ascertain whether the pixel was correctly assigned to a class or misassigned to
another class in a confusion matrix given in Table 2.
The errors are stated as commission and omission errors. Commission error
results from incorrect identification of a pixel, while omission error occurs when
we simply don’t recognize a pixel that we should have identified as belonging to a
particular class. The Quickbird true color image was used for visual reference.
The accuracy assessment of the two classifications shows that using the pro-
posed methodology has led to a significant increase in classification accuracies.
The accuracy achieved for the maximum likelihood classification is 67 and 79% in
urban and peri-urban areas, respectively.
The textural characteristics of the trees play a significant role in identification.
The smaller area of tree and the similar spectral signatures of grass, play grounds,
and lawns make it imperative to integrate textural parameters in classification
schemes. In order to achieve this, a selective smoothing procedure was adopted
on high resolution imagery to isolate and characterize individual trees. In the
Laplacian maxima filtering, this pre-processing step is crucial since this requires
 Remote Sensing

Table 2 Accuracy of the image classifications

Multi-spectral classification Laplacian maxima filtering

error
error

Palm
error
error

Palm
Urban Urban

accuracy
Mapping
accuracy
Mapping

Omission
Omission

Non-palm
Non-palm
area area

Commission
Commission

Palm 68 32 32% 33% 51.12% Palm 97 3 3% 5% 92.38%

Non-palm 33 67 33% 32% 50.75% Non-palm 5 95 5% 3% 92.23%
Total 101 99 67.5% Total 102 98 96%

error
error
error
error

Palm
Palm

Peri-urban

accuracy
accuracy

Mapping
Mapping

Omission
Omission

Non-palm
Non-palm

Commission
Commission

area
Palm 82 18 18% 24% 66.12% Peri-urban Palm 99 1 1% 3% 96.12%
Non-palm 24 76 24% 18% 64.40% area Non-palm 3 97 3% 1% 96.04%
Total 106 94 79% Total 102 98 98%
105
106 D. Al Ajmi and Saif ud din

a processing system to encapsulate the image content. Feature extraction results

are directly related to the performance of the initial enhancement and smoothing
stage. The methodology helped in accurate classification of date palm trees in
urban and agriculture areas. An accuracy of 96% is achieved for the urban area
and 98% for the agriculture area. The advantages of this methodology are clearly
the higher accuracies for similar scenes and simple processing. The problems
related to shadowing and irregular strands are solved to a larger extent, using this
approach. The image processing was carried out using PCI Geomatica software.
However, it is possible to utilize any other software for image processing (Figs.
43–48).
It was demonstrated that the proposed methodology surpasses the conventional
maximum likelihood classification in terms of palm tree identification in the study
area. The methodology used is precise in palm tree classification in both the urban
and agricultural areas. The classification in this instance is simple and straightforward,
and the subjective nature of data interpretation which relies to a great deal on the
interpreter’s perception is greatly minimized. These high accuracies in classifica-
tion are also due to the fact that there is no tree undergrowth and the uniformity of
palm tree crown sizes.
The analyses show that the palm trees are highlighted as blobs and the crown
patterns are very clearly segregated from other vegetation types in the area. The
Laplacian blob maxima coincide with the center of the tree. The success of the blob
detection technique in spatial mapping of palm trees in urban areas of Kuwait
shows potential of the methodology to be used in other Middle Eastern countries
with similar climatic and vegetation patterns.

Fig. 43 Wuickbird image, Kuwait city

Remote Sensing 107

Non Palm
Trees

Non Palm
Trees

Palm Trees

Fig. 44 Maximum likelihood classified image, Kuwait city

Fig. 45 Date palm tree blobs using laplaccian maxima (inverted LUT), Kuwait city

7.4 Application IV – Subsidence in Oil Fields

Consequent to the 1991 Gulf War there was extensive damage to the oil facilities.
After the war large-scale development was carried out leading to further expansion
of exploration and production facilities. The remote sensing technique has been
108 D. Al Ajmi and Saif ud din

Fig. 46 Quickbird image showing farm in north Kuwait

Non Palm Trees

Palm Trees

Non Palm
Trees

Fig. 47 Maximum likelihood classified image, showing farm in north Kuwait

successively employed in oil field areas to map the micro-elevation changes. Radar
interferometry for elevation mapping has been successfully used. The interferomet-
ric techniques give very precise measurements using highly correlated radar images
[39]. Surface sagging and subsidence in Oil fields have been reported from all over
the world [73–80]. The subsidence is usually gradual and often so small that it
remains undetected in a conventional ground survey. The surface subsidence in oil
Remote Sensing 109

Fig. 48 Date palm tree blobs using laplaccian maxima (inverted LUT), showing farm in north
Kuwait

fields has caused enormous financial losses due to disruption of production and loss
of infrastructure. Studies have been carried out to evaluate the impact of subsidence
in oil fields in California, Netherlands, Ekofisk – Norway, Brent – United Kingdom
[81]. Shell-operated production in the areas has been badly affected due to surface
subsidence. The InSAR technique is used in Kuwait to map the elevation changes
within the oil fields. Subsidence is typically explained by tectonic activity, but in
areas like Burgan Oil field, Kuwait, where there is little evidence of tectonic activ-
ity in the historical past, it is believed to be caused by mass balance adjustments.
The mass balance can lead to sagging of the surface and slow subsidence as hydro-
carbon reserves are removed from the host formation. The host formation under-
goes compaction due to over burden and loss of pressure from within the formation.
This can also reactivate sub-surface geological structures, which may compromise
the integrity of the reservoir seal, resulting in natural migration of hydrocarbons to
other formations. Furthermore, the compaction consequent to subsidence may lead
to reduced porosity which can lower the production levels in a formation.
Since the launch of the ERS-1 satellite by the European Space Agency (ESA) in
1991, the topic of Synthetic Aperture Radar (SAR) interferometric processing of
signals has gained attention in the remote sensing community all over the world.
The interferometric technique using SAR data is a very precise measurement tech-
nique to measure the deformation effect of the earth’s surface with sub-centimeter
accuracy—it is reported to be accurate even at the millimeter level [73, 80, 82]. The
radar interferometry is an alternative to the conventional stereoscopic method for
extracting topographic and deformation information. Synthetic aperture radar is an
110 D. Al Ajmi and Saif ud din

imaging radar device which images the radar backscatter of the earth’s surface over
large areas with a spatial resolution of 10–20m. InSAR uses the change in phase for
mapping the relative position of a point in space. In radar interferometry the change
in phase is measured between the backscatter in two different SAR images of the
same area taken from slightly different positions or at different times. The second
set of SAR data is taken with the same perspective and is co-registered and com-
bined into an interferogram. Pixel-wise phase difference is measured which indi-
cates the change in relative position of a pixel.
The basic idea of radar interferometry is based on the interference of electro-
magnetic waves. The interferometric technique relies on processing of two SAR
images of the same area from slightly displaced passes at different times or from
two antennas on the same platform at the same time. In the present communication
we are using images from slightly displaced passes at different times. The SAR is
an active remote sensing system where the sensor acts as transmitter and receiver
antenna during the image acquisitions.
Six SAR scenes corresponding to 09.05.1996, 12.05.1996, 13.05.1996,
29.05.1996, 10.03.1999, and 29.3.1999 from ERS 1 and 2 have been analyzed to
map the elevation changes. The data of 29.3.1999 was not helpful since it had a
perpendicular baseline of −1,504m which leads to loss of coherence. The other five
data scenes were used. The 1996 scenes are a tandem pair and were most useful to
give the topographic phase information, since the deformation, atmospheric and
noise phases are eliminated due to their insignificant contribution in the phase infor-
mation. The topographic information from the SAR data showed very good correla-
tion with the DEM data derived from the ground information. The phase information
is used in interferometry which gives millimeter-level vertical precision [83].
Because of such a fine vertical resolution it had been utilized in mapping subsid-
ence, tectonic deformation before, during and after earthquakes, inflation and defla-
tion of volcanoes due to magma movements [73]. It utilizes the correlation of
small-scale roughness to provide a coherent scattered signal, because of its high
sensitivity, care is taken so that the signals are not decorrelated due to the long time
interval of passes, vegetation growth and anthropogenic activity. In interferometry
the relative positional accuracy of each pixel is very high.
The coherence value of the three images in the Burgan oil field area shows high
coherence values of about 0.8–0.9. High coherence indicates low phase noise and
can give very precise estimates, while in lower values of coherence the phase noise
is pronounced and possibilities of precise deformation estimation are less. Figure 49
shows three image pairs showing different coherence values and subsidence infor-
mation. All the three show good coherence in the oil field area. The coherence
values are low in the western part of the scene probably due to the presence of a
sand sheet. In the present study repeat pass radar interferometry was used and inter-
ferograms are formed by two image signals, I1 and I2, which can be mathematically
expressed as:

{I = I I *
1 2 }
= A1e ij1 A2 e − ij2 = A1 A2 e i(j1 −j2 ) = A e iΦ ,
Remote Sensing 111

29.05 29.05
N N
W E W E
29 S 29 S

28.95 28.95
Legend
Legend 0.0 mm
28.9 1.0 28.9
Burgan Oil Field Burgan Oil Field
28.85 28.85

28.8 28.8

−50.0 mm
28.75 0.0 28.75
47.8 47.85 47.9 47.95 48 48.05 48.1 47.8 47.85 47.9 47.95 48 48.05 48.1
a b
The correlation (a) and subsidence (b) using scene dated 12.5.1996 and 13.5.1996
29.05 29.05
N N
W E W E
29 S 29 S

28.95 28.95
Legend
Legend 0.0 mm
28.9 1.0 28.9
Burgan Oil Field
Burgan Oil Field
28.85 28.85

28.8 28.8

28.75 0.0 28.75 −50.0 mm

47.8 47.85 47.9 47.95 48 48.05 48.1 47.8 47.85 47.9 47.95 48 48.05 48.1
a b
The correlation (a) and subsidence (b) using scene dated 12.5.1996 and 29.5.1996
29.05 29.05
N N
W E W E
29 S 29 S

28.95 28.95
Legend
Legend 0.0 mm
28.9 1.0 28.9
Burgan Oil Field Burgan Oil Field
28.85 28.85

28.8 28.8

0.0 28.75 −50.0 mm

28.75
47.8 47.85 47.9 47.95 48 48.05 48.1 47.8 47.85 47.9 47.95 48 48.05 48.1
a b
The correlation (a) and subsidence (b) using scene dated 12.5.1996 and 10.3.1999

Fig. 49 Correlation and subsidence images from Burgan Oil field

where * denotes the conjugate complex; A=A1A2 is interferogram amplitude;

F= 1– 2 is the interferogram phase.
For identical imaging geometry and surface I1=I2=0. But this is not observed in
practice, thus the phase difference is the sum of contributions from different proc-
esses [84] [85], which are expressed as:
Φ = Φtopo + Φdef + Φatm + Φn + Φdis ,
112 D. Al Ajmi and Saif ud din

where Φtopo is the phase contribution due to viewing of topography from different
perspectives, Fdef phase contribution due to movement along the line of sight, Fdis
phase contribution due to displacement, Fatm phase contribution due to atmospheric
effect, Fn is phase noise.
In order to determine the surface movement the assumption is taken that Forbit,
Fatm, Fdef, and Fn are negligible. Thus, Ftopo can be separated from Fdis. The separa-
tion of the topographic phase was done to yield the displacement phase. The inter-
ferometric phase component resulting from surface displacement is calculated
using expression:
Φdis = Φunwrappedphase − Φtopo .

The overall coherence values of the study area are in the range of 0.2–0.9. Burgan
Oil field shows a high coherence, which is good for subsidence measurements. The
tandem pair (12th and 13th May 1996) shows almost no subsidence in the highly
correlated area (Burgan Oil field). The data for 12th May 1996 was combined with
29th May data and no detectable subsidence was observed either. In the third case
the data of 12th May 1996 and 10th March 1999 were combined to yield maximum
displacement of 27mm in Burgan oil field area. The ASTER Mosaic of the study
area with corresponding subsidence rates shows a general overview of where sub-
sidence exists (Fig. 50).
The subsidence rates in the areas with low correlation are high, which may be
attributed to the sand movement, anthropogenic activities or simply are errors due
to atmospheric artifacts or baseline approximations. But in the area of Burgan oil
field where the correlation value is near unity, the displacement values seem to be
very reliable. The technique shows subsidence exists in the oil fields of Kuwait. The
diffused seismic activity in the region can be partly attributed to the mass balance
changes in the oil reservoirs, however the seismic activity is seen more often in the
Dahar area, where too it can be attributed to mass balance changes in the carbon-
ates. But in the case of Dahar and Managish oil field areas these changes are mani-
fested on the surface as sinkholes due to brittle failure of carbonate lithology. In the

a 29.05
N b 29.05
N
W E
W E
29 S 29 S

28.95 28.95
Legend
Burgan Oil Field 0.0 mm
28.9 28.9

28.85 28.85

28.8 28.8

28.75 −50.0 mm
28.75
47.8 47.85 47.9 47.95 48 48.05 48.1 47.8 47.85 47.9 47.95 48 48.05 48.1

Fig. 50 ASTER image of study area (left) and corresponding subsidence map (right)
Remote Sensing 113

case of Burgan the lithology is primarily clastics with intercalated shales. The
Burgan Formation, which is Albian in age, is 330m thick and forms part of the
Wasia Group [86, 87] which is overlain by Aruma, Hasa, and Kuwait Groups.
Burgan is the largest oil reservoir of Kuwait and it has produced for over seven
decades. The wells are under artesian conditions, which implies that the pressure in
the Formation is maintained. This pressure in host formations can be due to super-
incumbent strata of overlying formations. The mass balance adjustment due to
hydrocarbon extraction, leads to sagging of the surface and slow subsidence, as the
host formation adjusts to displacement of the hydrocarbon reserves. The infinitesi-
mal deformation over a period of years has apparently kept the oil facilities undis-
turbed and unaffected.

7.5 Application V – Sustainable Development of Fresh Water

Resources

The arid lands in the Middle-East have scarcity of water, which affects the human
population directly and can be a serious issue of conflict. The freshwater resources
are limited mostly to fossil water from the humid past; few shallow aquifers receive
recharge through precipitation in the arid ecosystem of the State of Kuwait and
adjacent countries, which calls for sustainable management of fresh water resources.
Kuwait has a two-tier aquifer system comprising the Dibdibah Formation of the
Kuwait Group at the top and the underlying Dammam Formation, which supports
the lower aquifer. The Kuwait Group has a shallow aquifer system of Quaternary
age, which is unconfined in nature. Lithologically it comprises silt and gravelly
sand. The Dammam Aquifer underlies the Kuwait Group, it is a chalky dolomitic
limestone of middle Miocene age, which is under confined conditions, however
upward leakages from the Dammam Formation are reported at places [88]. The
groundwater level varies from 90 m above mean sea level in the southwest to zero
at the Arabian Gulf in the northeast. The flow follows the regional northeastern dip.
The groundwater quality from these aquifers varies from brackish in the southwest
to highly saline in the northeast. The brackish water in the south and central part
with TDS around 4,000 ppm is used for irrigation and landscaping, but in the north
the water is highly saline with TDS exceeding 130,000 ppm [89, 90], which renders
it unfit for direct use. Only a few shallow aquifers receive groundwater recharge in
the middle-east, Raudatain Basin in Kuwait is one of them. The shallow aquifers
are irreparably damaged by over exploitation in the arid lands.
Quantification of hydrological budget is extremely difficult over large spatial and
temporal domains through direct observations, therefore remote sensing technology
has been extensively used to estimate various critical factors like landuse – land-
cover, runoff, evaporation, evapotranspiration [91–99]. In the present study remote
sensing technology has been used to estimate recharge to fresh water shallow aqui-
fers for their sustainable management in arid ecosystems. The fresh water in north-
ern Kuwait occurs as a lens in the Raudatain area. This lens floats over the highly
114 D. Al Ajmi and Saif ud din

saline groundwater in the northern part of the country. The Raudatain fresh water
lens is used for bottling drinking water, however the recharge to the fresh water lens
has not been estimated earlier. The area lies to the north of Kuwait city (Fig. 51).
Geologically the area is upper Miocene–Lower Pleistocene age and belongs to
the upper Dibdibah Formation of the Kuwait Group. Lithologically the formation
consists of coarse grained pebbly sand with thin intercalations of clayey sand and
clay, pebbles, cobbles, gravel, and conglomerate. The assemblage is similar to a fan
deposit. The surface geology shows consolidated calcritic deposit, which is helpful
in combating aeolian and fluvial erosion to a large extent, but impedes recharge to
the shallow aquifers (Fig. 52).
The Dammam Formation underlies the Kuwait Group. It is Middle Miocene in
age, lithologically comprising inter-bedded marine marls, limestone, and clays. The
thickness of this formation varies from 30 to 101m. The limestone horizons are
karstified in the Dammam Formation. The yield from the karstified limestone hori-
zons of the Dammam Formation estimated in the adjacent Al Hasa area in eastern
Saudi Arabia vary between 12 and 14 m3s−1 [100]. The study area shows signs of sec-
ondary salinity along the drainage channels and low infiltration of bedrocks (20 cmh−1)
[101] as is evident from seasonal playas that are formed after a rain event. The
entire State of Kuwait is dissected by a network of channels, most of which are
confined to the Dibdibah Formation (Kuwait Group) of Upper Miocene–Lower
Pleistocene age [101]. The drainage mapping of northern Kuwait demarcates the

Iraq Iran

Kuwait

Qatar 46.6 46.8 47 47.2 47.4 47.6 47.8 48 48.2 48.4 48.6
30.2 30.2
UAE
Saudi Arabia N
an

30 30
m

W E
O

n
me 29.8 S 29.8
Ye

29.6 29.6

29.4 29.4
Kuwait
Kuwait City
29.2 29.2

29 29

28.8 28.8

28.6 Study Area 28.6

28.4 28.4
46.5 47 47.5 48 48.5

Fig. 51 Location map

Remote Sensing 115

46.6 46.8 47 47.2 47.4 47.6 47.8 48 48.2 48.4 48.6

30.2
N

30
W E

S 29.8

29.6

29.4

29.2

LEGEND 29
Aeolian Sand
Alluvium
Desert Floor Deposits
28.8
Sabkah Deposits
Upper Dibdibah Formation
Lower Dibdibah Formation
Upper Member of Fars Formation 28.6
Lower Member of Fars Formation
Undifferentiated Fars and Ghar Formations
28.4
0 50 Kilometers

Fig. 52 Simplified geological map (After Saif uddin, et al., 2007)

catchment boundary of the Raudatain Basin using Landsat ETM data. The
Raudatain Basin shows centripetal drainage with low gradient. The hydrological
conditions suggest that the Raudatain lens received some recharge over the years,
which could be a fraction of actual drop from sky, possibly due to the development
of secondary salinity and playa over the recent past which has reduced the perme-
ability of the Dibdiba Formation to a greater extent (Fig. 53).
The recharge of the basin has been estimated by integration of the precipitation
data with the geology, geomorphology, and hydrological parameters [44].
The precipitation was estimated through interpretation of Tropical Rainfall
Measuring Mission (TRMM) data and scenes where rainfall events are observed
over the study area in the year 2003 have been selected. The rainfall rates in mil-
limeters were computed for each month over the basin (Table 3).
The DEM is created from SRTM data. The relief setting of the watershed is an
indicator of runoff/recharge potential. The terrain slope is used to estimate the
transmission loss, since channel slope affects the depth and duration of inundation.
116 D. Al Ajmi and Saif ud din

Fig. 53 Drainage map of Raudatain and Umm al Aish, north Kuwait

Table 3 Precipitation data over the study area

Accumulated rainfall in millimeters
Latitude Longitude January February March April May November December
29.5 47 18.1726 10.5877 21.9088 5.6055 1.2054 20.3955 21.2853
29.5 47.25 20.5399 7.8053 22.7921 5.7284 1.2279 25.3884 23.9906
29.5 47.5 27.3494 11.7285 24.2056 7.3434 1.8269 25.4238 25.6062
29.5 47.75 25.0254 2.6466 27.2732 6.9157 2.5216 20.1604 25.0226
29.5 48 30.6996 5.4489 23.9708 6.8856 0.8461 11.5332 28.8292
29.75 47 20.6627 7.2271 20.3377 7.1416 0.5861 20.4274 29.1998
29.75 47.25 25.8573 12.1315 24.653 7.048 0.8313 17.8435 24.4987
29.75 47.5 30.8797 11.0816 25.5302 6.7578 0.957 24.5263 27.944
29.75 47.75 36.1809 8.7045 27.9115 6.7853 1.5195 21.7922 26.2875
29.75 48 30.0956 11.6422 26.587 6.3967 1.7616 17.3794 29.0265
30 47 25.6618 8.4404 21.9767 8.7315 0.3186 23.5925 29.5724
30 47.25 27.3636 9.0023 22.8833 8.5124 0.5973 19.6751 25.3681
30 47.5 32.3377 7.0743 27.0348 9.7716 0.4416 21.4028 27.9787
30 47.75 33.1183 5.5189 26.0576 8.9359 0.1746 11.111 28.8616
30 48 31.3326 6.9656 27.3083 9.9083 0.4941 17.3641 28.9998
Mean monthly 27.68514 8.40036 24.69537 7.497847 1.02064 19.86771 26.8314
Net yearly rainfall 115.9985
Remote Sensing 117

W E
30
S
40 Elevation in Meter AMSL
240
29.8 220
200
180
160
29.6 140
120
100
80
60
29.4 40
20
0
47 47.2 47.4 47.6 47.8 48

Raudatain Basin 0 20 kilometer

00.20.4

Fig. 54 Digital Terrain Model of Raudatain Basin, north Kuwait

The slope also reflects the bed size material; the central depression in Raudatain
Basin shows a gentle slope and finer materials (playa) (Fig. 54).
The Natural Resources Conservation Service (NRCS) method of the United
States Department of Agriculture [102] was used to calculate the initial losses in
the sub-basin. The actual estimation of recharge is a complex process which takes
into account the initial losses and transmission losses. Initial losses occur in the
sub-basin before runoff reaches the stream networks, whereas the transmission
losses begin when the water is channeled through the stream network. The initial
losses are largely related to infiltration, surface soil type, landuse activities, eva-
potranspiration, interception and surface depression storage. Whereas the trans-
mission losses are primarily due to infiltration throughout the stream network.
The transmission losses are related to the channel geometry, upstream flow vol-
ume, duration of flow, bed size material, sediment load, and temperature. In this
study the transmission losses in the basin were computed using a regression
model developed by Walters [103] for similar climatic and lithological conditions
in Saudi Arabia. Walters in his model has modified the antecedent moisture index
equation developed by Linsley et al. [104] and in this model the soil ability to
accept infiltration was taken rather than a measure of soil moisture. Another
assumption inherent in the modification is that each flood would restore the soil
moisture to unity. The precipitation values are taken as cumulative monthly val-
ues calculated from TRMM data. The morphometric parameters of the micro
watershed are calculated. To have a synoptic view and to integrate the thematic
information that is required for recharge estimation in the micro watershed the
thematic files on topography, geology, geomorphology, and precipitation were
generated.
118 D. Al Ajmi and Saif ud din

The Landsat scene of path row 165–039 of March 2003 was taken since the tonal
contrast was good between the drainage channels, bare sand, lithology due to
higher moisture saturation. The cloud coverage was less than 5%.
The precipitation data calculated from TRMM over the basin were used to esti-
mate the recharge in the basin. The actual recharge to the underlying aquifer was
considerably less, since it equals the precipitation minus initial and transmission
losses. The initial losses occur in the basin before the runoff reaches the drainage
network, while there are transmission losses during the channeling of the water
through the drainage channels. The initial losses are largely related to the infiltra-
tion, surface soil type, landuse/landcover, evaporation, evapotranspiration, and
interception. The NCRS method used to estimate the initial losses suggest that
runoff in the basin will occur after rainfall exceeds an initial abstraction (Ia) value.
Rainfall excess, Q in the NRCS method is related to the effective precipitation
(P–Ia), through a maximum potential retention value, S. Thus, Q can be expressed
by equation, where all the measurements are in inches:

( P − I a )2
Q= .
(P − Ia + S )

The maximum potential retention S is a function of an empirical curve number

(CN) and is expressed as equation:

⎛ 1,000 ⎞
S=⎜ − 10⎟ .
⎝ CN ⎠

The initial abstraction as suggested by Gheith and Sultan (2002) for similar condi-
tions is 20% of maximum potential retention.
The CN is a function of antecedent moisture condition (AMC), the landuse, the
hydrologic condition, and soil type. In the study area the antecedent moisture con-
dition is believed to be low, since the rain events are rare and usually months apart.
The landuse and hydrologic conditions were classed as Natural Desert landscaping
and desert shrub coverage less than 10% of ground. Most of the basin is covered by
upper Dibdiba Formation and gravelly desert floor deposits which are localized to
the paleo channels with infiltration rates of 15 and 35 cm h−1, respectively [101].
The playa deposit at the center of the depression shows extremely low infiltration
of 4 cm h−1.
The soil type distribution in the area is related to the geomorphology and geol-
ogy of the area. Thus, we considered the curve number (CN) of the Dibdiba
Formation, desert floor deposit and playa in the basin. Here the sum of respective
CN values each weighted in proportion to its area is used to compute a composite
CN. The composite CN is 64 in the Raudatain Basin. Considering the aridity of the
area, extremely bare surface and high hydraulic conductivity of the desert floor
deposits, which are localized in the paleo drainage channels, it is assumed that
recharge to the alluvial aquifers can be approximated by transmission losses from
Remote Sensing 119

infiltration. These transmission losses are controlled by channel geometry, upstream

flow, duration of flow, lithology, sediment load, and temperature [105]. Evaporation
losses are inconspicuous given the short duration of rain event and the cloudy con-
ditions that prevail throughout the rain event.
The regression model [103] was used in the present case to estimate the trans-
mission losses. This model was developed for southwestern Saudi Arabia and for
identical climatic conditions and landform, thus we have used it to estimate the
transmission losses.
V1 = 0.1.3VA0.872 ,

where V1 =transmission loss in first mile; VA =upstream flow volume (acre-feet).

The transmission loss was calculated for each mile keeping the VA as remaining
volume from last mile.
The annual precipitation in millimeters over the entire watershed was converted
into the net volume by multiplying it with the basin area, which is estimated at
333.9596 MCM for the year 2003. The initial losses are estimated using the NCRS
method. The accumulated average annual initial loss is 78.43% of the net precipita-
tion, while the transmission losses are calculated by the Walters regression model
as 14.86% of the net precipitation. The net estimated annual recharge to the shallow
aquifer/lenses is estimated as 22.38 MCM. The monthly estimated water balance in
the study area is given in Table 4.
The Raudatain lens is tapped by 25 wells with a cumulative capacity of 4.5 mil-
lion liters per day. But due to an increase in TDS from 550 to 1,150 mgl−1 during
1964–1972 [106] the production of these wells was considerably reduced and pres-
ently it is used for bottling water. The lowering of the abstraction rate has resulted
in lowering of TDS levels of Raudatain water to about 600mgl−1 [107], which sug-
gest that some component of natural recharge does exist in the area. The 22.38
MCM of estimated recharge and the portion of 49.6485 MCM of transmission loss
annually in Raudatain Basin can be used to provide artificial recharge through
recharge wells to shallow aquifers and lenses. Construction of the recharge well
within the channels can reduce the transmission loss of 49.6485 MCM and a sig-
nificant portion of this component can be recharged to the shallow fresh water lens.
At certain point locations within the Kuwait Group the aquifer shows correlation

Table 4 Monthly water balance estimates (in MCM)

January February March April May November December Annual
Drop from 79.70552 24.18464 71.09798 21.5863 2.938423 57.19913 77.2476 333.9596
sky
Initial losses 59.77914 22.49171 54.38995 20.50699 2.938423 44.04333 57.78121 261.9307
Transmission 13.15141 1.692925 11.69562 1.079315 0 9.437856 12.59136 49.64848
losses
Monthly 6.774969 0 5.012408 0 0 3.717943 6.875036 22.38036
recharge
120 D. Al Ajmi and Saif ud din

between salinity distribution and channel ordering, i.e., low salinity had been
observed in higher order channels in a down stream direction, which suggest that
fresh recharge has taken place in the area. It is proposed to construct recharge wells
on higher order channels which will facilitate additional recharge to the shallow
fresh water lens, resulting in reduced transmission and evaporation losses, which
take place due to formation of playas on the surface.
The shallow aquifer and fresh water lens in Raudatain can be used to supply
fresh water for domestic consumption and can partly reduce pressure on desalina-
tion activity which is increasing the salinity at the coast in the Arabian Gulf. This
study projects the possibility of sustainable water resource development in the
desert state of Kuwait.

References

1. Curran PJ (1985) Principles of remote sensing. Longman Publishing Co., London, pp

208–224
2. Jensen JR (1986) Introductory digital image processing: a remote sensing perspective.
Printice-Hall, Englewood Cliffs, NJ, pp 1–9, 177–229
3. Lillesand TM, Kiefer RW (1987) Remote sensing and image interpretation, 2nd ed. Wiley,
New York, p 721
4. Pratt WK (1992) Digital image processing, 2nd ed. Wiley, New York, p 698
5. Sabins, Floyd F (1987) Remote sensing – principles and interpretation. Freeman and
Company, New York, p 449
6. Lee CC (1983) Elimination of redundant operations for fast sobel operator. IEEE Trans Syst
Man Cybern, vol SMC-13, No 3, pp 242–245
7. Pizer SM (1987) Adaptive histogram equalization and its variations. Comput Vision Graph
Image Process 39:355–368
8. Chavez PS Jr, Berlin GL, Sower LB (1982) Statistical method for selecting LANDSAT MSS
ratios. J Appl Photograph Eng 8:23–30
9. Roberts LG (1965) Machine perception of three-dimensional solids in optical and electro-optical
information processing. JT Tippett et al. (ed) MIT Press, Cambridge, MA, pp 159–197
10. Kirsch R (1971) Computer determination of the constituent structure of biomedical images.
Comput Biomed Res 4(3):315–328
11. Lewis TO, Odell PL (1971) Estimation in linear models. Prentice-Hall, Englewood Cliff, NJ
12. Prewitt JMS (1970) In: Lopkoin BS, Rosenfeld A (eds) Object enhancement and extraction in
picture processing and psychopictorics. Academic Press, New York
13. Press WH, Flamery BP, Tenkolsky SA, Vetterling WT (1988) Numerical recipes in C: The art
of scientific computing. Cambridge University Press, New York
14. Jahne B (1995) Digital image processing, concept, algorithms and scientific applications, 3rd
ed. Springer-Verlag, Berlin, p 383
15. Rosenfeld A, Kak AC (1982) Digital image processing, 2nd ed., vol. I and II. Academic Press,
Orlando
16. Short NM (1982) The LANDSAT tutorial workbook. Basics of satellite remote sensing.
National Aeronautics and Space Administration, Reference Publication 1078, Washington, DC
17. Loeve M (1955) Probability theory. D. Van Nostrand Company, Princeton, NJ
18. Castleman KR (1977) Digital image processing. Prentice Hall, Englewood Cliff, NJ
19. Kwarteng AY, Al-Ajmi D (1997) Satellite remote sensing applications in the State of Kuwait.
Kuwait Institute for Scientific Research, Kuwait, p 101
Remote Sensing 121

20. Saif ud din, Al Dousari A, Literathy P (2008) Evidence of hydrocarbon contamination from
the Burgan Oil Field, Kuwait – interpretations from thermal remote sensing data. J Environ
Manage 86:605–615
21. Saif ud din, Al Dousari A, Ramdan A, Al Ghadban AN (2008b) Site specific precipitation
estimate from TRMM data using bi-linear weighted interpolation technique – an example
from Kuwait. J Arid Environ. doi:10.1016/[Link].2007.12.013
22. McClain EP, Pichel WG, Walton CC (1985) Comparative performance of AVHRR-based
multichannel sea surface temperatures. J Geophys Res 90:11587–11601
23. Richter R (1990) A fast atmospheric correction algorithm applied to Landsat TM. Int J
Remote Sens 11:159–166
24. Pearce AF, Prata AJ, Manning CR (1992) Comparison of NOAA/AVHRR-2 sea surface tem-
peratures with surface measurements in coastal waters. Int J Remote Sens 10:37–52
25. Qin Z, Karnieli A, Berliner P (2001) A mono-window algorithm for retrieving land surface
temperature from Landsat TM data and its application to the Israel–Egypt border region. Int
J Remote Sens 22(18):3719–3746
26. Liang S (2001) An optimization algorithm for separating land surface temperature and emis-
sivity from multispectral thermal infrared imagery. IEEE Trans Geosci Remote Sens
39(2):264–274
27. Sobrino JA, Li ZL, Stoll MP, Becker F (1996) Multi-channel and multi angle algorithm for estimat-
ing sea and land surface temperature with ASTER data. Int J Remote Sens 17(11):2089–2114
28. Gillespie AR, Rokugawa S, Matsunaga T, Cothern J, Hook S, Khale A (1998) A temperature
and emissivity separation algorithm for advanced spaceborne thermal emission and reflection
radiometer (ASTER) images. IEEE Trans Geosci Remote Sens 36:1113–1126
29. Jimenez-Munoz JC, Sobrino JA (2003) A generalized single channel method for retrieving
land surface temperature from remote sensing data. J Geophys Res 108.
doi:10.1029/2003JD003480
30. Sobrino JA, Jimenez-Munoz JC, Paolini L (2004) Land surface temperature retrieval from
Landsat TM 5. Remote Sens Environ 90:434–440
31. Markham BL, Barker JL (1986) Landsat MSS and TM post calibration dynamic ranges,
atmospheric reflectance and at satellite temperature. EOSAT Landsat Technical Notes 1
(Lanham, Maryland: Earth Observation Satellite Company), pp 3–8
32. Schneider K, Mauser W (1996) Processing and accuracy of landsat TM data for lake surface
temperature measurement. Int J Remote Sens 16:2111–2124
33. Hurtado E, Vidal A, Caselles V (1996) Comparison of two atmospheric correction methods
for Landsat TM thermal band. Int J Remote Sens 17:237–247
34. Literathy P (1992) Environmental consequences of the Gulf War in Kuwait: impact on water
resources. Water Sci Technol 26:21–30
35. Litherathy P (1993) Considerations for assessment of environmental consequences of the
1991 Gulf War. Marine Pollut Bullet 27:349–356
36. Kwarteng AY (1998) Multitemporal remote sensing data analysis of Kuwait’s oil lakes.
Environ Int 24(1/2):121–137
37. Kwarteng AY (1999) Remote sensing assessment of oil lakes and oil-polluted surfaces at the
Greater Burgan oil field, Kuwait. Int J Appl Earth Obs Geoinf 1(1):36–47
38. Al-Dousari A (2001) Analysis of change in the tarcrete layer on the desert surface of Kuwait
using satellite imagery and cell-based modeling. PhD thesis remote sensing & GIS. Boston
University (unpublished)
39. Saif uddin, Al-Dousari A, Al-Ghadban AN, Aritoshi M (2006) Use of interferometric tech-
niques for detecting subsidence in the oil fields of Kuwait using Synthetic Aperture Radar
Data. J Petrol Sci Eng 50:1:1–10
40. Al-Senafy MN, Viswanathan MN, Senay Y, Sumait A (1997) Soil contamination from oil
lakes in northern Kuwait. J Soil Contaminat 6(5):481–494
41. Robertson FR, Fitzjarrald DE, Kummerow CD (2003) Effect of uncertainty in TRMM precipi-
tation radar path integrated attenuation on interannual variations of tropical oceanic rainfall.
Geophy Res Lett 30(4). doi:10.1029/2002GL016416
122 D. Al Ajmi and Saif ud din

42. Anon (2005) Tropical rainfall mission office – introduction. [Link]/lba/eng/trmm/

[Link]
43. Bowman KP (2005) Comparison of TRMM precipitation retrievals with rain gauge data from
ocean buoys. J Climate 18:178–190
44. Saif ud din, Al Dousari A, Al Ghadban AN (2007) Sustainable freshwater resources manage-
ment in northern Kuwait – A remote sensing view from Raudatain basin. Int J. Appl. Earth
Obs. Geoinfo 9:21–31
45. Brito AE, Chan SH, Cabrera SD (2003) SAR image superresolution via 2-D adaptive extrapo-
lation. Multidimen Syst Signal Process 14:83–104
46. Kuhel H, Sacchi MD (2003) Least squares wave equation migration for AVP/AVA inversion.
Geophysics 68:262–273
47. Gribbon KT, Bailey DG (2004) A novel approach to real-time bi-linear interpolation. Second
IEEE international workshop on electronic design, test and applications, p 126
48. Foley TA (1987) Weighted bicubic spline interpolation to rapidly varying data. ACM Trans
Graph 6(1):1–18
49. Arnold DN, Boffi D, Falk RS (2002) Approximation by quadrilateral finite element. Math
Comput 71:909–922
50. Nespor V, Sevruk B (1999) Estimation of wind induced error of rainfall gauge measurements
using a numerical simulation. J Atmos Oceanic Technol 16:450–464
51. Habib E, Krajewski WF, Nespor V (1999). Numerical simulation studies of rain gauge data
correction due to wind effect. J Geophy Res 104:19,723–19,734
52. Krajewski WF, Ciach GJ, McCollum JR, Bacotiu C (2000) Initial validation of Global pre-
cipitation climatology project monthly rainfall over United States. J Appl Meteorol
39:1071–1086
53. Serra YL, A’Hearn P, Freitag HP, McPhaden MJ (2001) ATLAS self-siphoning rain gauge
error estimates. J Atmos Oceanic Technol 18:1989–2002
54. Bell TL, Kundu PK (2003) Comparing satellite rainfall estimates with rain gauge data: Optimal
strategies suggested by spectral model. J Geophy Res 108 (D3) 4121, pp ACL 7-1- 15
55. Sawaya WN (1986) Dates of Saudi Arabia. Regional agriculture and water research center.
Ministry of Agriculture and Water, Riyadh, Saudi Arabia, p 200
56. Zang Y (2001) Texture – integrated classification of urban treed areas in high resolution color
infrared imagery. Photogram Eng Remote Sens 67(12):1359–1365
57. Barnsley MJ, Barr SL (1996) Inferring urban land use from satellite sensor images using
kernel-based spatial reclassification. Photogram Eng Remote Sens 62:949–958
58. Gao J, Skillcorn D (1998) Capability of SPOT XS data in producing detailed land cover maps
at the urban-rural periphery. Int J Remote Sens 19:2877–2891
59. EEC (1997) Decision No. C (97)2690. Decision of commission of the financing of an action
concerning counting of the number of Olive trees in member states.
60. Howard JA (1991) Remote sensing of forest resources. Chapman & Hall, London
61. Marceau DJ, Gratton DJ, Fournier RA, Fortin JP (1994) Remote sensing and the measurement
of geographical entities in a forested environment. 2. The optimal spatial resolution. Remote
Sens Environ 49(2):105–117
62. Woodcock CE, Strahler AH (1987) The Factor of scale in remote sensing. Remote Sens
Environ 21:311–332
63. Dreyer P (1993) Classification of land cover using optimized neural nets on spot data.
Photogramm Eng Remote Sens 59:617–621
64. Anys H, Bannari A, He DC, Morin D (1994) Texture analysis for the mapping of urban areas
using airborne MEIS-II images. Proceedings of the first international airborne remote sensing
conference and exhibition, Strasbourg, France, vol 3, pp 231–245
65. St-Onge BA, Cavayasm F (1997) Automated forest structure mapping from high resolution
imagery based on directional semivariogram estimates. Remote Sens Environ 61:82–95
66. Culvenor DS, Coops N, Preston R, Tolhurst K (1999) A spatial clustering approach to auto-
mated tree crown delineation. Proceedings of the international forum on automated interpreta-
tion of high spatial resolution digital imagery for forestry, 10–12 February, Victoria, BC,
Canada, Natural Resources Canada, pp 67–80
Remote Sensing 123

67. Coops N, Culvenor DS (2000) Utilizing local variance of simulated high spatial resolution
imagery to predict spatial pattern of forest strands. Remote Sens Environ 71:248–260
68. Wulder M, Niemann KO, Goodenough DG (2000) Local maxima filtering for the extraction
of tree locations and Basdal area from high spatial resolution imagery. Remote Sens Environ
73:103–114
69. Granlund GH (1999) The complexity of vision. Signal Process 74:101–126
70. Kay S, Leo O, Peedel S, Giardino G (2000) Computer assisted recognition of olive tree in
digital imagery. Space Applications Institute JRC of the European Commission, Ispra, Italy
71. Karantzalos K (2003) Combining anisotropic diffusion and alternating sequential filtering for
satellite image enhancement and smoothing. Proceedings international conference of image
and signal processing for remote sensing IX, SPIE, Barcelona, Spain
72. Alvarez L, Lions PL, Morel JM (1992) Image selective smoothing and edge detection by
nonlinear diffusion. SIAM – JNA 29:845–866
73. Biegert EK, Berry JL, Oakley SD (1997) Oil field subsidence monitoring using spaceborne
interferometric SAR – A Belridge 4-D Case History. Proceeding annual meeting of the
American association of petroleum geologists, Dallas, April’97
74. Carne C, Massonnet D, King C (1996) Two examples of the use of SAR interferometry on
displacement fields of small spatial extent. Geophys Res Lett 23(24):3579
75. Briole P, Massonnet D, Delacourt C (1997) Posteruptive deformation associated with the 1986–
1987 and 1989 lava flows of Etna detected by radar interferometry. Geophys Res Lett 24(1):37
76. Fielding EJ, Blom RG, Goldstein RM (1998) Rapid subsidence over oil fields measured by
SAR interferometry. Geophys Res Lett 25:3215
77. Galloway DL, Hudnut KW, Ingebristen SE, Philips SP, Peltzer G, Rogez F, Rosen PA (1998)
InSAR detection of aquifer system compaction and land subsidence, Antelope Valley, Mojave
Desert, California: Water Resources Research, vol 34, pp 2573–2585
78. Galloway DL, Jones DR, Ingebristen SE (1999) Land subsidence in the United States: US
Geological Survey Circular 1182, pp 1–117
79. Galloway DL, Jones DR, Ingebristen SE (2000) Measuring land subsidence from space; US
Geological Survey fact sheet 051–00, pp 1–4
80. Hoffmann J, Zebker HA, Galloway DL, Amelung F (2001) Seasonal subsidence and rebound
in Las Vegas valley, Nevada, observed by synthetic aperture radar interferometry. Water
Resour Res 37(6):1551–1566
81. Tuttle M, Ehrismann J, Hulshof B. (1998) Detection and monitoring of surface subsidence
using synthetic aperture radar interferometry in Yibal Oil Field, Sultanate of Oman. http://
[Link]/papers/patzek/[Link]
82. Gray AL, Farris-Manning PJ (1993) Repeat pass interferometry with airborne synthetic aper-
ture radar. IEEE Trans Geosci Remote Sens 31(1):189–191
83. Hoffmann J (2003) The application of satellite radar interferometry to the study of land sub-
sidence over developed aquifer systems. PhD dissertation, Department of Geophysics,
Stanford University, USA, pp 1–211
84. Zebker HA, Werner CL, Rosen PA, Scott H (1994) Accuracy of topographic maps derived
from ERS – 1 interferometric radar. IEEE Trans Geosci Remote Sens 32(4):823–836
85. Ferretti A, Prati C, Rocca F (2000) Non-linear subsidence rate estimation using permanent
scatter in differential SAR Interferometry. IEEE Trans Geosci Remote Sens
35(5):2202–2212
86. Khalaf FI, Gharib IM, Al-Hashash MZ (1984) Types and characteristics of the recent deposits
of Kuwait, Arabian Gulf. J Arid Environ 7:9–33
87. El-Baz F, Al-Sarawi M, Al-Shamlan AA (2000) Atlas for the State of Kuwait from Satellite
Images. Kuwait Foundation for the Advancement of Sciences, pp 1–145
88. Abusada SM (1981) Interpretation of environmental isotope data of Kuwait groundwater.
Proceeding first symposium on the future development of water resources in the Gulf and
Arabian Peninsula. J Arab Gulf Arab Peninsula Stud 2:8–92
89. Burdon DJ, Al Sharhan A (1968) The problem of the paleokarstic Dammam Limestone
Aquifer in Kuwait. J Hydrol 6:358–404
 124 D. Al Ajmi and Saif ud din

90. Himida IH, El-Yaqubi SA (1979) Hydrogeological and hydrogeochemical aspects of the
mail groundwater fields in Kuwait, Arabian Gulf first geological congress of the middle east,
Ankara, Turkey, pp 1–4
91. Anderson JR, Hardy EE, Roach JT, Witmer RE (1976) A landuse and landcover classifica-
tion system for use with remote sensor data. USGA Professional Paper 964, p. 28
92. Beven KJ, Wood EF, Sivapalan M (1988) On hydrological heterogeneity: Catchment mor-
phology and catchment response. J Hydrol 100:353–375
93. Menenti M, Bastiaanssen WGM, van Eick D, El Karim MAA (1989) Linear relationships
between surface reflectance and temperature and their application to map actual evaporation
of groundwater. Adv Space Res 9(1):165–176
94. Menenti M (1993) Understanding land surface evapotranspiration with satellite multispectral
measurements. Adv Space Res 13(5):89–100
95. Bastiaanssen WGM, Hoekman DH, Rorbrling RA (1994) A methodology for assessment of
surface resistance and soil water storage variability at mesoscale based on remote sensing
measurements. IAHS Special Publications 2, Institute of Hydrology, Wallingford, UK, p 66
96. Bastiaanssen WGM, Menenti M, Dolman AJ, Feddes RA, Pelgrum H (1996) Remote
sensing parameterization of meso scale land surface evaporation. In Raschke E (ed)
Radiation and water in the climate system: remote measurements. Springer-Verlag,
Berlin, pp 401–429
97. Bastiaanssen WGM, Menenti M, Feddes RA, Holtslag AAM (1998). A remote sensing sur-
face energy balance algorithm for land (SEBAL). Part 1: formulation. J Hydrol
212–213:198–212
98. Bastiaanssen WGM, Pelgrum H, Wang J, Ma Y, Moreno JF, Roerink GJ, van der Wal T
(1998). A remote sensing surface energy balance algorithm for land (SEBAL).: Part 2:
Validation J Hydrol 212–213:213–229
99. Bastiaanssen WGM, Noordman EJM, Pelgrun H, Davids G, Thoreson BP, Allen RG (2005)
SEBAL Model with remotely sensed data to improve water resource management under
actual field condition. J Irrigat Drain 131(1):85–93
100. Edgell HS (1997) Aquifers of Saudi Arabia and their geological framework. Arab J Sci Engg
22(1C):4–31
101. Al Sulaimi J, Khalaf FJ, Mukhopadhyay A (1997) Geomorphological analysis of paleo
drainage systems and their environmental implications in the desert of Kuwait. Environ Geol
29(1–2):94–111
102. SCS (1985) National engineering handbook, Section 4: Hydrology, US Department of
Agriculture, Soil Conservation Service, Engineering Division, Washington DC
103. Walters MO (1990) Transmission losses in arid regions. J Hydraul Eng ASCE 116:129–138
104. Linsley RK, Kohler MA, Paulhus JLH (1975). Hydrology for engineers. McGraw Hill, New
York, NY
105. Gheith H, Sultan M (2002) Construction of a hydrologic model for estimating Wadi runoff
and groundwater recharge in the Eastern Desert, Egypt. J Hydrol 262:36–55
106. Senay Y (1977) Groundwater resources and artificial recharge in Raudatain waterfield.
Ministry of Electricity and Water, Kuwait, p 35 (unpublished)
107. Kwarteng AY, Viswanathan MN, Al Senafy AN, Rashid T (2000) Formation of fresh ground-
water lenses in northern Kuwait. J Arid Environ 46:137–155

View publication stats