GEOSTATIONARY SATELLITE AND RAINGAGE NETWORK DATA

FOR RAINFALL ESTIMATION

By

KE-SHENG CHENG

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL

OF THE UNIVERSITY OF FLORIDA IN

PARTIAL FULFILLMENT OF THE REQUIREMENTS

FOR THE DEGREE OF DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1989

ACKNOWLEDGMENTS

I wish to express my sincere appreciation to the

chairman of my supervisory committee, Dr. S.F. Shih, for

his guidance and assistance throughout my graduate program.

I also thank Dr. E.P. Lincoln, Dr. A.R. Overman, Dr. M.A.

Uman, and Dr. M.C. Yang, members of the supervisory

committee, for their help and encouragement.

Special thanks are extended to Dr. J.R. Carr of the

Department of Geological Science, University of Nevada-Reno

and Dr. D.E. Myers of the Mathematics Department,

University of Arizona for providing valuable research

papers and the Co-Kriging program. Thanks are also given

to Mr. O. Lanni and Mr. C. Bryant for their help with

various aspects of this research.

Thanks are also given to South Florida Water Manage-

ment District for providing rainfall data and other

information.

I thank my parents for their endless love,

expectation, and firm support. I also wish to express my

deepest appreciation to my wife, Yen-Ching, for

accompanying me through all the time of happiness and

frustration.

TABLE OF CONTENTS

Page

ACKNOWLEDGMENTS . . . . . . . . .

LIST OF TABLES . . . . . . . . .

LIST OF FIGURES . . . . . . . . .

ABSTRACT . . . . . . . . . . .

CHAPTERS

1 INTRODUCTION . . . . . . . .

. ii

S vi

S viii

. xi

The Influence of Spatial Rainfall Variation

on Hydrological Modeling . . . . .

Satellite Data for Rainfall Estimation . .

The Use of Co-Kriging . . . . . .

Objectives . . . . . . . . .

2 REVIEW OF LITERATURE . . . . . . .

Life-history Methods . . . . . .

Stout, Martin, and Sikdar Method . . .

Griffith-Woodley Method . . . . .

The life history technique . . . .

The streamlined technique . . . .

Rainfall apportionment technique . . .

Scofield-Oliver Method . . . . . .

Pattern Classification Methods . . . .

Wu, Weiman, and Chin Method. . . . .

Weiss-Smith Method.. . . . . . .

Kriging/Co-Kriging Approach of Estimation. .

3 CHARACTERISTICS OF GOES SATELLITE SYSTEM

AND DATA . . . . . . . .

GOES Instruments and Data Characteristics

VISSR Instrument Characteristics . . .

4 STUDY AREA AND DATA SETS . . . . .

Hourly Rainfall Data . . . . . .

GOES Infrared Data . . . . . .

Pixel Registration of GOES Satellite Data.

S. 29

S. 29

S. 30

S. 34

iii

. .

. .

. .

5 THEORY OF REGIONALIZED VARIABLES

--KRIGING/CO-KRIGING . . . . .

Second-Order Stationarity and Intrinsic

Stationarity ... . . . . . .

Computing the Experimental Variogram . .

Commonly Used Variogram Models and

Their Properties . . . . . .

The Continuity . . . . . . .

The Range of Influence . . . . .

The anisotropies . . . . . .

Parameter Inference of Selected Variogram

Model . . . . . . .

Effect of Oversampling of GOES IR Data

on Variogram Estimation . . . .

Optimum Estimation by Kriging.. . . .

Optimum Estimation by Co-Kriging . . .

6 RAINFALL ESTIMATION--CO-KRIGING APPROACH .

Phase I--Pixels Classification . . .

Delineating Cloud-covered Area . . .

Determining Classification Features. .

Mathematical Concepts of Pattern

Classification . . . . . .

Classification of Cloud Pixels Using

Bayes Optimal Classifier . . . .

Phase II--Estimation of PXLRN for Ungaged

Pixels Using Co-Kriging Method . . .

Formulation of Co-Kriging Equation

for PXLRN Estimation . . . . .

Random Field Analysis . . . . .

Estimation of variograms . . . .

Dimensionless variograms . . . .

Dimensionless cross-variograms .

Criteria for Model Verification and

Results Evaluation . . . . .

Cross-Validation . . . . .

. 45

. 46

. 46

. 47

. 50

. 50

. 50

. 51

. 51

. 52

. 54

. 57

. 57

. 57

. 58

. 61

. 63

. 75

75

78

79

79

102

*. 109

S. 112

7 RESULTS OF RAINFALL ESTIMATIONS. . .

Verification of the Model Validity and

Results Evaluation . . . . .

Possible Sources of Error. . . .

Incomplete Observation of Rainfall .

Pixel Registration . . . . .

Calculation of PXLRN of Gaged Pixels .

Estimation of Dimensionless Variograms

Comparing Estimates of Test and Training

Areas . . . . . . . .

* . 113

* . 113

* . 126

S . 126

. . 126

. . 127

. . 127

. . 128

8 CONCLUSIONS .. . . . . . .. .

. 132

. .

APPENDICES

A GROUND-TRUTH PXLRN, CO-KRIGING ESTIMATES, AND

VARIANCES OF ESTIMATION ERRORS . . .. 135

B EFFECTS OF THE RANDOM FIELD VARIANCE

ESTIMATION ON THE CO-KRIGING ESTIMATES AND

ESTIMATION ERROR VARIANCES . . . .. 148

C EXPLANATION OF GOES IR DATA RESOLUTION . . 151

D EXPERIMENTAL DIMENSIONLESS VARIOGRAMS OF

PXLRN, CCTT, AND (PXLRN+CCTT). . . ... 152

REFERENCES . . . . . . . . ... ... 165

BIOGRAPHICAL SKETCH. . . . . . . . .. 169

LIST OF TABLES

Table Page

2-1 Rainfall rate as a function of echo growth

trend for south Florida echoes . . ... 12

2-2 Radiance and textural features used by

Wu et al. (1985) . . . . . ... 25

4-1 Historic storm events . . . . ... 36

4-2 Geographic coordinates of raingages inside

the study area . . . . . . ... .37

5-1 Commonly used variograms . . . ... 48

6-1 Misclassification rates of pattern classi-

fication technique--all images . . ... 65

6-2 Expanding/contracting status of GOES

satellite images . . . . . ... .68

6-3 Statistics of PXLRN of expanding and

contracting images . . . . . ... .70

6-4 Misclassification rates of pattern classi-

fication technique--expanding images . .. 72

6-5 Misclassification rates of pattern classi-

fication technique--contracting images. . 73

6-6 Estimated parameters of dimensionless

variograms of different storms . . ... 99

7-1 Verification of assumptions and estimations

for Storm-79215 . . . . . ... 114

7-2 Verification of assumptions and estimations

for Storm-82252 . . . . . ... 115

7-3 Verification of assumptions and estimations

for Storm-82264 . . . . . .... 116

7-4 Verification of assumptions and estimations

for Storm-82265 . . . . ...... 117

7-5 Verification of assumptions and estimations

for Storm-82269 . . . . . ... 118

7-6 Regressions of estimated PXLRN on

ground-truths . . . . . . ... 124

7-7 Comparison of ratios of absolute error in

test and training areas . . . ... 130

7-8 Comparison of variances of estimation

errors in test and training areas ... . 131

vii

LIST OF FIGURES

Figure

2-1 An idealized schematic history of a cloud

on a satellite image and its associated

echo . . . . . . . . . .

2-2 Cloud and echo area relationships for GOES

infrared data over Florida . . . .

2-3 Digital enhancement curve for GOES infrared

images . . . . . . . . .

2-4 Schematic diagram of Scofield-Oliver method

for rainfall estimation . . . .

3-1 Oversampling of GOES infrared data . .

4-1 The study area and raingage locations .

4-2 Life history of Storm-79215 . . . .

4-3 Life history of Storm-82252 . . . .

4-4 Life history of Storm-82264 . ....

4-5 Life history of Storm-82265 . . . .

4-6 Life history of Storm-82269 . . . .

5-1 Commonly used variogram models . . .

6-1 A pixel group used to determine values of

classification features . . . . .

6-2 Evolution of average STDDEV . . .

6-3 Histogram of PXLRN of expanding and

contracting images . . . . . .

6-4 A cloud system showing precipitation-free

and precipitating areas . . . . .

6-5 Experimental variograms of PXLRN of

Storm-82269 .. . . . . . . .

Page

. 9

S 10

19

. 21

. 32

. 35

. 40

S. 41

S. 42

S. 43

S. 44

S. 49

S. 59

. 66

. 69

S. 76

. 80

viii

6-6 Changes of the exponential variogram when

w is independent of t . . . . ... 82

6-7 Changes of the exponential variogram when

a is independent of t . . . . . .. 83

6-8 Changes of the exponential variogram when

both w and a are functions of t . . .. 84

6-9 Average of PXLRN versus standard deviation

of PXLRN . . . . . . . . 85

6-10 Dimensionless variograms of

Storm--79215 . . . .

6-11 Dimensionless variograms of

Storm--82252 . . .

6-12 Dimensionless variograms of

Storm--82264 . ...

6-13 Dimensionless variograms of

Storm--82265 . . .

6-14 Dimensionless variograms of

Storm--82269 . . .

6-15 Dimensionless variograms of

Storm--79215 . . .

6-16 Dimensionless variograms of

Storm--82252 . . . .

6-17 Dimensionless variograms of

Storm--82264 . .

6-18 Dimensionless variograms of

Storm--82265 . ...

6-19 Dimensionless variograms of

Storm--82269

PXLRN for

PXLRN for

PXLRN for

PXLRN for

PXLRN for

CCTT for

CCTT for

CCTT for

CCTT for

CCTT for

* . 89

. . 90

. . 91

. . 92

. . 93

. . 94

S. 95

* . 96

. . 97

. . . . . . . . 98

6-20 Fitted dimensionless variograms of PXLRN

6-21 Fitted dimensionless variograms of CCTT

6-22 Dimensionless variograms of standardized

(PXLRN+CCTT) for Storm--79215 . . .

6-23 Dimensionless variograms of standardized

(PXLRN+CCTT) for Storm--82252 . . .

6-24 Dimensionless variograms of standardized

(PXLRN+CCTT) for Storm--82264 . . .

ix

S. 100

. 101

S. 104

S. 105

. 106

6-25 Dimensionless variograms of standardized

(PXLRN+CCTT) for Storm--82265 . . .

6-26 Dimensionless variograms of standardized

(PXLRN+CCTT) for Storm--82269 . . .

6-27 Fitted dimensionless variograms of

standardized (PXLRN+CCTT) . . . .

6-28 Dimensionless cross-variograms . . .

7-1 Ground-truth PXLRN versus Co-Kriging

estimates (Storm-79215) . . . . .

7-2 Ground-truth PXLRN versus Co-Kriging

estimates (Storm-82252) . . . . .

7-3 Ground-truth PXLRN versus Co-Kriging

estimates (Storm-82264) . . . . .

7-4 Ground-truth PXLRN versus Co-Kriging

estimates (Storm-82265) . . ....

7-5 Ground-truth PXLRN versus Co-Kriging

estimates (Storm-82269) . . . . .

7-6 Ground-truth PXLRN versus Co-Kriging

estimates (all storms) . . . . .

S. 107

S. 108

S. 110

. 111

S. 119

S. 120

. 121

S. 122

S. 123

S. 125

Abstract of Dissertation Presented to the Graduate School

of the University of Florida in Partial Fulfillment of the

Requirements for the Degree of Doctor of Philosophy

GEOSTATIONARY SATELLITE AND RAINGAGE NETWORK DATA

FOR RAINFALL ESTIMATION

By

KE-SHENG CHENG

December 1989

Chairman: S.F. Shih

Major Department: Agricultural Engineering

Geostationary Operational Environmental Satellite

(GOES) infrared images and rainfall measurements from

raingage network were used together to estimate pixel-

averaged rainfall (PXLRN) for ungaged locations. Five

storm events that occurred over Central Florida during

August 1979 and September 1982 were studied.

A threshold value of satellite-derived cloud-top

temperature was used to identify cloud-covered area. GOES

satellite images were divided into two groups based on the

expanding/contracting status of cloud systems over the

study area. Clouds of expanding stage were found to have

higher cloud-top temperature variation over space than do

clouds of contracting stage. Most heavy rainfalls

(rainrate > 0.2 inch/hr) occurred during the expanding

stage. A statistical pattern classification technique

xi

using Bayes optimal classifier was applied to assign cloud-

covered pixels to either a precipitation-free pixel or

precipitating pixel. Two classification features, coldest

cloud-top temperature (CCTT) and standard deviation

(STDDEV) of cloud-top temperature within a group of nine

pixels, were used for classification. Different cutoff

rainrates were used to define precipitation-free and

precipitating pixels. The study area includes a training

area with high raingage density and a test area with low

raingage density. Using cutoff rainrate of 0.02 inch/hr,

more than 70% of pixels in the training area were correctly

classified for both expanding and contracting images, and

59% and 64% of pixels in the test area were correctly

classified for expanding and contracting images,

respectively.

Co-Kriging method of estimation was used to estimate

PXLRN for ungaged precipitating pixels using PXLRN and CCTT

of gaged pixels. In addition to PXLRN estimation, Co-

Kriging method also provides variance of estimation error.

The unbiasness of Co-Kriging was well satisfied. Paired t-

test showed no significant difference between PXLRN

estimates and their ground-truths. Correlation

coefficients of PXLRN estimates and their ground-truths

vary from 0.62 to 0.85 for different storms. Overall

correlation coefficient is 0.86.

xii

CHAPTER 1

INTRODUCTION

Rainfall-runoff relationship is one of the most

important topics in hydrologic research. Rainfall

information is the input of many hydrologic models.

However, a raingage network only provides scattered

rainfall measurements and in general is not capable of

describing the spatial variation of rainfall needed for

distributed hydrological models. Recently, many

hydrologists have been investigating the application of

satellite data to rainfall estimation. One of their major

concerns is how the weather satellite data can be used to

help describe the spatial distribution of rainfall so that

rainfall, along with other hydrologic parameters (for

example, soil moisture content and permeability) and

spatial data (such as land-use information), can be

included in a geographic database. Then a geographic

information system (GIS) can be used to help develop

distributed hydrological models.

The Influence of Spatial Rainfall Variation

on Hydrological Modeling

The success of flood routing heavily depends on the

information of rainfall distribution over the watershed.

2

Conventional methods of hydrologic flood routing require

the information of average rainfall intensities over the

watershed. Rainfall data measured by raingages are usually

used to obtain the areal average rainfall. However, the

spatial variation of rainfall has been found to have a

remarkable influence on the behavior of the runoff

hydrograph (Wilson et al., 1979). It was also established

that, at least for drainage basins with areas less than

about 3,000 km2, the spatial and temporal rainfall pattern

is important for the determination of runoff, and that the

climatological raingage spacing is inadequate for the

definition of hydrologically useful spatial precipitation

patterns (Georgakakos and Kavvas, 1987; Hamlin, 1983;

Nicks, 1982). Collier (1985) pointed out that, even in

cases where the total depth of precipitation is correctly

estimated by an existing raingage network, and even when

the temporal character of the precipitation is accurately

recorded, there could be serious errors in the total

volume, peak, and time-to-peak of the estimated hydrograph

when the spatial pattern of precipitation was not properly

represented in the input data. O'Connell and Clark (1981)

suggested that it might be more beneficial to seek a better

representation of the spatial variation in rainfall and its

effect on streamflow response, and an improvement in the

structure of real-time forecasting models, than to expand

effort in solving runoff estimation problems.

3

Satellite Data for Rainfall Estimation

Rainfall estimation utilizing remotely sensed data has

drawn the attention of many hydrologists and meteorologists

since the early 1960s. Over the last three decades several

rainfall estimation techniques utilizing weather satellite

data have been developed. In areas where data from more

than one sensor (e.g., raingage, radar, and satellite) are

available, it appears that use of all the data should

produce better precipitation estimates than those obtained

with the data from each sensor (Georgakakos and Kavvas,

1987). This appears so because various sensors have

different capabilities of spatial coverage and measure

different variables that are related to precipitation. The

raingage measures precipitation depth at ground level

directly. The satellite such as Geostationary Operational

Environmental Satellite (GOES) measures irradiances of the

cloud tops. Compared with data from raingage network,

satellite data not only have larger spatial coverage but

also have much higher spatial resolution. Recent research

has shown that data from raingage network is insufficient

for the estimation of a time-space stochastic model of

rainfield and additional remotely sensed data are necessary

(Karr, 1986).

The Use of Co-Kriging

The spatial variations of rainfall and extracted

variables of satellite data and the relationships between

4

rainfall and these variables are far too complex to be

described by usual analytical functions. Although one

would usually expect the colder cloud pixels to have more

rainfall than other warmer pixels, Whitney and Herman

(1979) found that satellite IR black-body temperatures

correlated poorly with point observations of hourly

rainfall. The cloud-top temperature alone is either not

enough for rainfall estimation or the relationship between

rainfall and cloud-top temperature is too erratic and, of

course, not monotonic.

In order to consider the spatial variations of both

raingage measurement and cloud-top temperature and their

spatial correlation, Co-Kriging, an estimation method of

regionalized variables, is thus employed for rainfall

estimation.

Objectives

The main objective of this study was to use both

raingage measurements and geostationary satellite data to

estimate rainfall for ungaged areas. The specific

objectives were

1. to study the spatial distribution characteristics

of raingage measurements;

2. to develop a technique for identifying expanding

and contracting clouds;

3. to identify precipitation-free and precipitating

areas on satellite images;

5

4. to develop dimensionless variograms for both

pixel-averaged rainfall and cloud-top temperature;

5. to develop dimensionless cross-variogram of

pixel-averaged rainfall and cloud-top temperature;

6. to develop a Co-Kriging approach using both

satellite-derived cloud-top temperatures and

raingage measurements to estimate rainfall for

ungaged locations; and

7. to compare the results of Co-Kriging estimates

with the raingage measurements.

CHAPTER 2

REVIEW OF LITERATURE

Since the launch of the first meteorological satellite

(TIROS I, 1 April 1960), hydrologists and meteorologists

have been looking for techniques that estimate rainfall

utilizing information provided by satellites. Various

investigations have shown that rainfall from convective

clouds can be identified and estimated objectively based on

cloud-top temperature fields derived from GOES infrared

images. This chapter summarizes some of these techniques.

They are classified into two categories based on the

characteristics of each technique: life-history and pattern

classification methods.

Rainfall data and satellite-derived cloud-top

temperatures are essentially spatial data. Recently, the

Kriging/Co-Kriging method has been applied to analyze many

spatially distributed hydrologic parameters including

rainfall. Some of these works are also briefly described

in this chapter.

Life-history Methods

Stout, Martin and Sikdar Method

From June through September 1974, a multination Global

Atmospheric Research Program Atlantic Tropical Experiment

6

7

(GATE) was conducted with two major objectives, as stated

by Kuettner (1974):

1. To estimate the effects of smaller-scale tropical

systems on large-scale circulations.

2. To advance the development of numerical modeling

and prediction methods.

Cloud information, including total cloud amount, low

cloud amount, height of base of lowest cloud, etc., was

collected from GATE shipboard radars. Ground measurements

of rainfall were also made by radars. During GATE the

Synchronous Meteorological Satellite (SMS-1) provided

visible (VIS) and infrared (IR) wavelength images at an

interval of 30 minutes or less. Stout, Martin and Sikdar

(1979) studied GATE data and proposed the following

equation for estimation of rainfall rate:

Rv = agAc + al(dAc/dt) (2-1)

where Rv is volumetric rainrate (13/t) for a particular

cloud, Ac is the area of the cloud at time t, dAc/dt is

cloud area change, and ag and al are empirical

coefficients.

In sequences of satellite pictures matched to

simultaneous radar images the areas and volumetric rainrate

of individual cumulonimbus clouds were measured throughout

the cloud lifetimes. The cloud area was approximated by a

single threshold contour of temperature for infrared data

or brightness for visible data. The thresholds were 200

W/m2 and -280C for visible and infrared data, respectively.

8

The coefficients a0 and al were calculated from the

combined measurements by least square regression of

satellite cloud area on radar volumetric rainrate. The

standard error of estimate was 60% of the volumetric

rainrate for visible, 75% for infrared (Stout et al.,

1979). Wylie (1979) presents evidence of a need to adjust

the coefficients if this scheme is used to estimate

rainfall in higher latitudes. For this method there is no

provision accounting for rainfall from clouds which are not

cumulonimbus, and each cumulonimbus cloud of a given size

and growth rate produces an identical rainfall.

Griffith-Woodlev Method

The fundamental premise of this method was that the

time behavior of clouds on the satellite images and the

corresponding radar echoes approximate the simple model as

shown in Figure 2-1 in which both clouds and radar echoes

grow to a maximum and then decay.

Griffith et al. (1978) studied complete histories of

several Florida cloud systems collected during GATE using

GOES satellite and radar observations. A threshold of

-200C was used to define cloud areas for GOES infrared

data. Woodley et al. (1980) and Griffith (1987) proposed

three techniques of rainfall estimation using satellite

data for different purposes.

The life history technique. Empirical relationships

between satellite-measured cloud areas and radar echo areas

were derived and shown in Figure 2-2. Cloud areas are

tl te

SATELLITE:

o D

Ac2

RADAR :

Ac3

Ac4 = An

Ac5

Ac6

V 0

Ae2

Ae3

Ae5

t; time index

Ac: satellite, cloud-covered area

Ae: radar echo area Am: maximum cloud-covered area

Figure 2-1. An idealized schematic history of a cloud on a satellite image

and its associated echo.(from Griffith et. al, 1978).

0.2

E

0.1 -

0i 0

0,0 0.2 0,4 0.6 0.8 1,0 0.8 0.6 0.4 0.2 0.0

Ac/Am

Figure 2-2. Cloud and echo area relationships for GOES infrared

data over Florida (from Griffith et al., 1978).

11

determined from series of satellite infrared data with

threshold of -200C and then associated radar echo areas can

be determined from Figure 2-2.

Radar echo areas are assumed related to volumetric

rainfall in the form of the following equation (Griffith,

1987):

Rv(tj) = I(tj)*Ae(tj)*At* E [ai(tj)*bi]*103

1

(2-2)

where tj indicates the image date and time; Rv is

satellite-inferred rain volume (m3); I is rainrate (mm/hr);

Ae is inferred echo area (km2) defined by the 1 mm/hr

rainrate; At is the time interval between successive images

(h); ai is the fractional coverage of the cloud by colder

temperatures; bi is an empirical weighting coefficient

(also a function of temperature); and i denotes different

temperature layers in a satellite image. Using digitized

Miami WSR-57 data, Griffith and Woodley derived a function

as shown in Table 2-1 for the determination of rainrate I

(Woodley et al., 1980; Griffith et al., 1980).

The values of the weighting coefficient b were

empirically derived and are a function of cloud-top

temperatures T (expressed as digital IR count DN),

according to

b = exp(0.02667+0.01547*DN)/11.1249

154 < DN < 176 or -200C > T > -310C

Table 2-1 Rainfall rate as a function of echo growth trend

for south Florida echoes (From Griffith et al.,

1980)

Echo Growth Trend Rainfall Rate (mm/hr)

Increasing Echo Area

0.00 < Ae/Aem < 0.25 13.3

0.25 < Ae/Aem < 0.50 17.3

0.50 < Ae/Aem < 0.75 21.1

0.75 < Ae/Aem < 1.00 23.8

Maximum Echo Area

Ae/Aem = 1.00 20.7

Decreasing Echo Area

1.00 > Ae/Aem > 0.75 21.1

0.75 > Ae/Aem > 0.50 16.7

0.50 > Ae/Aem k 0.25 11.9

0.25 > Ae/Aem > 0.00 8.2

Ae = echo area defined by the 1 mm/hr rain rate;

Aem = maximum area an echo attains in its life cycle.

13

b = exp(0.11537+0.01494*DN)/11.1249

176< DN < 255 or -310C > T > -1100C

(2-3)

The life history technique of Griffith-Woodley method

has been widely tested in various convective situations and

different regions. The tests include thunderstorms in

Florida, Montana, and Venezuela; hurricanes in the western

Atlantic; and the thunderstorms and cloud clusters of GATE

from West Africa across the Atlantic Ocean to South America

(Barrett and Martin, 1981). In Florida over a large area of

105 km2 with 49 cases of hourly rainfall, the correlation

coefficient of SMS-1 infrared image rainfall with radar

measurements of rainfall was 0.65. The infrared data showed

a tendency to overestimate rainfall (Barrett and Martine,

1981). Shih (1989a, 1989b) used the life history technique

to estimate daily rainfall for six cloud systems which

occurred in Florida. Hourly rainfalls were estimated for

individual satellite images and then accumulated to obtain

daily rainfall. Estimated daily rainfalls were found well

correlated with raingage measurements (r2=0.82).

The streamlined technique. A major disadvantage of

the life history technique is that it is not configured to

make estimates of rainfall in real time. The cloud system

must be traced through its life cycle in order to determine

the maximum cloud area which is essential in the rain

calculation. The streamlined technique was designed to

overcome this disadvantage. It makes estimation only from

14

single images so that rain estimates can be made in real

time.

The streamlined technique assumes that the dependence

of the inferred rainfall on the system's stage in its life

cycle is small compared with the dependence on cloud area

alone. In the streamlined technique the following

identifications are made:

I = 16.7 mm/hr (2-4)

Am = Ac (2-5)

where Ac is the cloud area on the current image, and

Ae/Am = 0.016 Ac < 2,000 km2

Ae/Am = 0.047 2,000 < Ac < 10,000 km2

Ae/Am = 0.067 Ac > 10,000 km2

For the streamlined technique,

0.016

Rv = [16.7*( 0.047 )*Ac*At*Eaibi]*103 (2-6)

0.067

Shih (1989a, 1989b) estimated daily rainfall for six cloud

systems using the streamlined technique. Results of

streamlined technique do not match the raingage

measurements as well as those from the life history

technique.

Rainfall apportionment technique. The life history

and streamlined techniques estimate volumetric rainfall of

cloud systems. Although areal-averaged rainfall can be

determined by dividing the inferred volumetric rainfall

with the cloud areas, we are still unable to calculate

rainfall depth at each pixel. The apportionment technique

15

was designed to apportion the total volume of rainfall to

pixels of a cloud system.

The apportionment technique places half of the

calculated rain volume into pixels that constitute the

coldest 10% of the cloud area; the remaining half of

rainfall is apportioned over pixels whose temperatures fall

into the next warmer 40% of the cloud area. Only half of

the cloud area that was used to compute the rain volume

will contain a rain rate, and the coldest cloud tops have

the most rain.

The following equation is used for calculating rain

depth of pixels:

Di = ------- 10-3 (2-7)

n

where Di is the rainfall depth (mm) of the ith pixel; Rv is

the inferred volumetric rainfall for this image; bi is the

weighting coefficient for the temperature of the ith pixel;

Eb is the sum of the weighting coefficients with the index

n which is the total number of pixels of either the coldest

10% of the cloud area or the next warmer 40% of the cloud

area, and g is the area of a pixel in km2.

At this point, Eqs. 2-3 and 2-7 are worthy of a

careful examination. Since the inferred volumetric

rainfall (Rv), the sum of the weighting coefficients (Eb),

and also the area of a pixel (g), are all determined and

fixed before performing the calculation of Di, the rainfall

depth is apparently a linear function of bi within the

16

coldest 10% or the next warmer 40% of the cloud area.

Again, the weighting coefficient b is an exponential

function of digital IR count, DN, as represented in

Eq. 2-3. For two pixels, say pixels i and j, with

different DN values, DNi and DNj, and both within either

the coldest 10% or the next warmer 40% of the cloud area,

the following relationship holds:

(DNi-DNj)

bi/bj = 1.0156 154 < DNi,DNj < 176

(DNi-DNj)

bi/bj = 1.0151 176 < DNi,DNj < 255

(2-8)

Consequently,

(DNi-DNj)

Di/Dj K1 1.0156 154 < DNi,DNj < 176

(DNi-DNj)

Di/Dj = K2 1.0151 176 < DNi,DNj < 255

(2-9)

where K1, K2 are constants.

Accordingly, rainfall depth (of the apportionment

technique) within either the coldest 10% or the next warmer

40% of the cloud area is a strictly monotonic function of

the digital IR count, DN, of pixels.

Nergi et al. (1983) examined the life history and the

apportionment techniques using a single day's GOES infrared

images over Florida. The magnitude and variability of

several terms in Eq. 2-2 were analyzed. Several of their

results and conclusions were drawn as follows:

17

1. Due to the definition of clouds by the threshold of

-200C, a simple model of cloud growth and decay as depicted

in Figure 2-1 did not exist for most clouds. Eighty-two

percent of the clouds they studied existed for one hour

(two images) or less.

2. The derived parameters of fractional echo area

(Ae/Am) and rainrate were notably uncorrelated with the

calculated rain volume.

3. The maximum cloud area term Am was highly

correlated with the calculated rain volume (r=0.95).

4. Correlation between half-hourly estimates of area-

averaged rainrate and gauge amounts was poor (r=0.03).

(However, the method was used for the calculation of area

averaged rainrate from gauge measurements is not clear.)

5. The life history technique is not generally

applicable to clouds defined by the -200C isotherm.

Scofield-Oliver Method

This method was originally developed for estimating

half-hourly convective rainfall by analyzing the changes in

two consecutive GOES satellite images. Unlike the rainfall

apportionment technique of the Griffith-Woodley method

which must estimate the volumetric rainfall of a cloud

system for a single image, the Scofield-Oliver method

estimates convective rainfall at interested locations while

not estimating the rain volume of the cloud system.

The estimation scheme of this method is based on the

18

following considerations (Scofield and Oliver, 1977):

1. Bright clouds in the visible imagery produce more

rainfall than darker clouds.

2. Bright clouds in the visible and clouds with cold

tops in the IR imagery which are expanding in areal

coverage (in early and mature stages of development)

produce more rainfall than those not expanding.

3. Decaying clouds produce little or no rainfall.

4. Clouds with cold tops in IR imagery produce more

rainfall than those with warm tops.

5. Clouds with cold tops that are becoming warmer

produce little or no rainfall.

6. Merging of cumulonimbus (Cb) clouds increases the

rainfall rate of the merging clouds.

7. Most of the significant rainfall occurs in the

upwind (at anvil level) portion of a convective system.

The highest and coldest clouds form where the thunderstorms

are most vigorous and the rain heaviest.

A digital enhancement curve, the Mb curve as shown in

Figure 2-3, was used to produce enhanced IR images. The

technique can be represented in the form of a decision

tree:

Step 1. Locate active portion of the thunderstorm system.

The first decision is to determine if the area of

concern is under the active upwind portion of the

convective cloud system by using enhanced IR and

VIS imagery. Clues for locating the active portion

COUNT VALUE (INPUT)

000

WHITE

>-

w

-J

BLACK _

134

56.8

330

88.4 42.5 -3.4 -74.5

31.3 5.8 -19.7 -59.2

305 279 253' 214

TEMPERATURE

Figure 2-3. Digital enhancement curve for GOES infrared

images (from Scofield, 1987).

255

204

0

051

C

Z

.-4

153

m

-4

-4

1000

- 165 oF

-109 C

164 OK

20

of a thunderstorm system are listed in step 1 of

Figure 2-4.

Step 2. Compute half-hourly rainfall estimate from the

factors described in step 2 of Figure 2-4.

Estimates of convective rainfall are made by

comparing changes in two consecutive pictures using

both IR and high resolution VIS images.

Step 3. Compute the total half-hourly convective rainfall

estimate. The total half-hourly convective

rainfall can be estimated by summing those

estimates computed in factors 1, 2, 3 and 4 and

then multiplied by the moisture correction factor.

The assignment of particular rainfall rates to

categories of cloud-top temperature, growth, position,

texture and relationship with other convective systems is

based on standard gauge measurements of convective rainfall

over the central United States for one summer season. As

illustrated in Figure 2-4, the Scofield-Oliver method only

estimates rainfall which occurred at the active portion of

the thunderstorm system. Therefore, the number of

estimates per cloud are quite variable and depend on the

number of heavy rainfall signatures analyzed within the

particular cloud. Although the Scofield-Oliver method is

not suitable for studying the spatial variation of

rainfall, it can be used for flash flood detection and

forecasting. Also, even Scofield (1987, p.1783) claimed

STEP 1

Rainfall is computed only for the active portion

of the thunderstorm system.

The following are clues for helping to make decision.

Area of IR temperature gradient at upwind end of

anvil for a thunderstorm system in strong vertical

wind shear.

Center of the anvil with a tight, uniform IR tempe-

rature gradient around entire anvil for a thunder-

storm system with no vertical wind and shear.

Area near and under an overshooting top.

Portion of anvil that is brighter and/or more

textured.

Half of anvil bounded by edge which moves least

(comparison of the last two IR or VIS pictures).

Area near 'upper level' (500 200 mb) upwind end of

anvil.

Area near low level inflow.

Area under a radar echo.

STEP 2

Half-hourly rainfall estimates in mm are computed

from the following factors

Factor 1

Cloud-top temperature and cloud growth factor.

Determine amount that the coldest cloud tops

increased within half-hour and select estimate

amount according to grey shade.

>1/3 Coldest

>2/3' <2/3* <1/3"Lat Area tops

C Lat Lat or same decrease warmer

Med Gray

Lt Gray

Dk Gray

Black

Rpt Gray

White

(-32 to -41)

(-41 to -52)

(-52 to -58)

(-58 to -62)

(-62 to -80)

(below -80)

6

13

19

25

25-50

50

4

8

10

15

15-25

25

2

4

5

8

8-15

15

Figure 2-4. Schematic diagram of Scofield-Oliver method for

rainfall estimation (Scofield, 1987).

Divergence aloft factor

Determine estimate according to gray shade when strong

divergence is apparent.

Med Gray

Lt Gray Dk Gray Black Rpt Gray

White

4 8 10 15 15-25 25

I

Factor 2

Overshooting Top Factor. Estimate an additional 8 mm

for colder tops in the area of the overshooting tops.

I

Factor 3

Thunderstorm or Convective Cloud Line Merger Factor.

Estimate an additional 12.7 mm for colder tops in the

area of the merger.

Factor 4

Saturated Environment Factor. Estimate additional

amount according to gray shade and time of persistence.

Med Gray Lt Gray Dk Gray Black Rpt Gray white

>1 hr 5 5 5 5 8 8

>2 hrs 10 10 10 10 13 13

Factor 5

Moisture Correction Factor =

Precipitable Water Relative Humidity

STEP 3

Total Half-Hourly Convective Rainfall Estimates (mm) =

sum of factors 1 through 4 and multiplied by moisture

correction factor.

Figure 2-4. Continued

23

that "Once a meteorologist acquires a chemistry with the

methodology estimates become reproducible, reliable, and

useful." One without profound background in meteoro-logy

may encounter difficulties when trying to locate the active

portion of a convective cloud system and when considering

other factors.

Pattern Classification Methods

Wu. Weinman and Chin Method

Infrared and visible satellite data obtained from

summer storms, tropical storms and cyclones in the Gulf of

Mexico coastal region and near south Florida were analyzed

for the development of this method (Wu et al., 1985).

A pattern classification algorithm was developed to

derive rainfall rates from single visible and infrared

images or infrared images only. This method determines

rainfall in 20x20 km grid elements on the basis of radiance

and texture features. Radiance features are global

measurements of the 20x20 km grid elements containing

information about the overall characteristic of the

radiances within the grid element. Texture features are

measurements that concern the spatial distribution of the

radiances.

This method classifies rainfall into three classes

that correspond to operational radar precipitation levels:

0--no rain (0 < R < 0.02 inch/hr); 1--light rain (0.02 < R

< 0.2 inch/hr); and 2--heavy rain (0.2 inch/hr < R).

24

Rainfall rates were inferred from radar data that were

obtained simultaneously with the satellite imagery.

Rainfall rates, R (mm/hr), were computed from the measured

radar reflectives, Z (mm6/m3), by the relationship Z =

55R1.6.

As shown in Table 2-2, 24 features from each of the

visible and infrared images were employed to determine the

various rainfall classes. The radiance and texture

features were regarded as components of a vector X that

describes each 20x20 km grid element.

A feature selection procedure which ranked the

features in the order of their ability to separate the

rainfall rate classes was applied to several training data

sets. Wu et al. (1985) found that the visible radiance

features K=[l], [3], and [4] in Table 2-2 were frequently

significant discriminating features for the separation of

rain (class 0) and no rain (class (1+2) ); however, visible

texture features K=[ll], [14], and [15] were also

frequently significant discriminating features. The

infrared radiance features K=[1], [3] and [4] and texture

features K=[10], [12] and [15] were frequently significant

in identifying precipitating regions. Separating rainfall

classes 1 and 2 was mostly achieved with features [1], [3]

and [4] in the visible and in the infrared images.

A detailed inspection of the contrast texture feature

showed that values of the contrast feature varied signifi-

cantly in different stages of convective activity of cloud

Table 2-2.

Radiance and textural features used by Wu et al.

(1985I

Feature

n-mmr III

Radiance

Mean grey level

Standard deviation of grey level

Maximum grey level

Minimum grey level

Maximum/minimum grey level ratio

Grey level range

Texture

[7] Edge strength per unit area

[8]-[11] Maximum within 4 directions

[12]-[15] Mean of 4 directions

[16] Edge strength per unit area

[17]-[20] Maximum within 4 directions

[21]-[24] Mean of 4 directions

. pxlsprto

-d = pixel separation

2RG--Rberts gradient

MEAN-mean of grey level differences

CON-contrast of grey level differences

ASM--angular second moment of grey level

ENT--entropy of grey level differences.

VIi

RG[7]2 2

MEAN[8],CON[9] 2

ASM[10],ENT[11]

MEAN[12],CON[13] 2

ASM[14],ENT[15]

RG[16] 4

MEAN[17],CON[18] 4

ASM[19],ENT[20]

MEAN[21],O1N[22] 4

ASMr231.ENTr241

differences

dl

S IR

4

4

~.c~YVw, i",Vn4=

L J f L J

Features

26

systems. Thus, contrast was selected as an initial

screening feature which divides training data sets into

three categories, i.e., high, medium, and low contrasts.

The maximum likelihood classifier was used as the first

decision rule to discriminate between the no-rain, 0, and

the rain (1+2) classes. The features used for this first

classification were visible K = [1], [3], [4], [7], [8],

[10], [11], [12], [14], [15], [19], [23], [24] and infrared

K = [1], [3], [4], [10], [14], [15] for high contrast cases

and visible K = [1], [3], [4], [14], [15], [19], [23],

[24] and infrared K = [1], [3], [4], [10], [14], [15] for

low contrast cases. After determining that each grid

element was either in the "non-raining" (0) or the

"raining" (1+2) classes, those grid elements in which rain

was falling were subsequently divided into classes 1 or 2.

The features used for this second classification were the

visible and infrared radiance mean, maximum and minimum

grey levels, K=[l], [3] and [4]. Of their results, case

studies that employed data from both the visible and

infrared sensors correctly identified rainfall classes 0

and (1+2) in about 90% of the cases and identification into

classes 1 and 2 was correct in about 70% of the cases

studied. Data derived only from infrared images yield

correct identification of 0 and (1+2) classes in 85% of the

cases and identification of cases 0, 1 and 2 was correct in

65% of the cases.

Weiss-Smith Method

Similar to Wu et al. (1985), Weiss and Smith (1987)

used a statistical pattern classification technique to

identify intense rainfall regions occurring at low

frequency from higher frequency light and no-rainfall

regions occurring within large convective cloud systems

over the United States High Plains. GOES infrared data at

full resolution (7 km at satellite subpoint) were used to

calculate a variety of parameters including coldest cloud-

top temperature (CCTT), average cloud-top temperature

(AVGCTT), and the difference of AVGCTTs of two consecutive

images (DEL-AVGCTT). Three different rainrate cutoffs:

0.04, 0.16, and 0.40 inch/hr were used. The overall

performance of the classification technique increased as

the rainrate threshold was raised from 0.04 to 0.16 inch/hr

and finally to 0.40 inch/hr. The parameters CCTT, AVGCTT

and DEL-AVGCTT were found to be the best discriminant

parameters for rainfall classification. Light rain events

were not easily distinguishable from no-rain events,

suggesting that light rainfall may occur too randomly

within the cloud to be adequately classified.

Kriging/Co-Kriqing Approach of Estimation

Since Matheron proposed the Kriging method in the

early 1960's, the primary applications of this method have

been in mineral resource investigations. However, its

successes have attracted the attention of researchers in

diverse fields.

Delhomme (1978) applied the Kriging method to the

calculation of the average rainfall depth over a given

region. The mean variogram of 13 storms was first

calculated. Then a hypothesis was made by assuming that

variograms of different storms are proportional to each

other based on their variances of rainfall depth. Finally,

the Kriging system for every storm was solved to obtain the

optimal weights to be ascribed to the raingages. In

forming the Kriging system, the mean value of the variogram

between Xi, the location of any one raingage, and a point

describing the domain, usually the watershed, must be

obtained in advance.

Cho (1985) addressed the issue of precipitation as a

dynamic process resulting from interactions at various

scales in the atmosphere vs. precipitation as a stochastic

process. He showed that small perturbations of the large-

scale temperature field can lead to significant

precipitation events at the mesoscale. Randomness is

introduced due to our inability to measure accurately, and

continuously in space and time precipitation controlling

variables such as temperature. He concluded by

recommending the study of the precipitation process by

coupling the dynamic evolution of the atmosphere to a

stochastic description of the precipitation-controlling

variables (e.g., temperature field).

CHAPTER 3

CHARACTERISTICS OF GOES SATELLITE SYSTEM AND DATA

The Geostationary Operational Environmental Satellite

(GOES) system includes the satellite itself; the National

Environmental Satellite, Data and Information Service

(NESDIS) facility at Wallops Island, Virginia (that

receives the direct transmission and generates the relay

"stretched VISSR" transmission); and the ground systems at

NESDIS and at the Space Science and Engineering Center

(SSEC) of the University of Wisconsin.

GOES Instruments and Data Characteristics

The original GOES instrument was the Visible and

Infrared Spin Scan Radiometer (VISSR), which was an

outgrowth of the spin scan radiometer flown aboard several

of the Applications Technology Satellite (ATS) series of

NASA research satellites. The GOES satellites are spin

stablized satellites operated at the geostationary altitude

of nominally 35,800 km (19,312 nml). The nominal satellite

spin rate is 100 revolutions per minute. VISSR is a two

channel instrument, containing a visible channel for high

resolution daytime cloud imaging and an IR water vapor

Inm: nautical mile

30

window (11 pm) to allow daytime and nighttime imaging and

the estimation of cloud height as a function of cloud-top

temperature. The next step in the development of GOES

instrumentation is the VISSR Atmospheric Sounder (VAS)

flown aboard GOES satellites from GOES-4 on. The VAS

instrument features a multispectral IR imaging capability

useful for deriving vertical profiles of atmospheric

parameters such as temperature and water vapor as well as

multispectral imagery.

VISSR Instrument Characteristics

The VISSR instrument consists of scanning system,

telescope, and Infrared and Visible sensors. The scanning

system consists of a mirror that is stepped mechanically to

provide North to South viewing, while the rotation of the

GOES satellite provides West to East scanning. Each step

of the mirror causes a change of 192 microradians (mr) in

the scan angle, representing a distance of 6.87 km, or

nominally 7 km, near nadir. A sequence of 1,821 scans is

performed to provide a "full disk" view from just beyond

the Northern Earth horizon to just beyond the Southern

Earth horizon. At the satellite rotation rate of 100 rpm,

18.21 minutes are required to complete one full North to

South view of the Earth.

The scanning mirror reflects the received radiation

into a sixteen-inch diameter telescope. A fiber-optics

bundle is used to couple the telescope to eight VIS

31

detectors (sensitive to the 0.54 to 0.70 pm band). The

fiber optics bundle is configured such that each of the

eight VIS sensors has a 25 (W-E) by 21 (N-S) mr

instantaneous field of view (IFOV). The IFOV provides a

ground resolution of 0.9 km, or nominally 0.5 nm. Thus the

system provides eight parallel VIS data lines per West to

East scan, covering the 7 km (4 nm) band scanned by each

step of the scanning mirror. The IFOV of the IR detectors

is 192 mr and thus the resolution near nadir is 6.87 km, or

about 4 nm. The output from VIS and IR detectors are

digitized on board the satellite and transmitted down to

Earth in real-time. The visible data is sampled every 2

ps, which yields visible samples spaced at increments of

satellite rotation of 20.9 mr, a near nadir spacing of 0.75

km. The IR data is sampled every 8 ps, which yields IR

samples spaced at increments of satellite rotation of 83.8

mr, a near nadir spacing of 3.0 km. Since the IR detector

IFOV is 192 mr, the IR data is therefore oversampled along

the scan. The oversampling of the IR data leads to the

designation of IR data as "4*2" IR data (4 nm resolution N-

S, 2 nm resolution W-E). Figure 3-1 illustrates the

oversampling of GOES IR data.

The real-time transmission of GOES VISSR data are

received at the NESDIS ground station and processed by the

Synchronous Data Buffer (SDB). The SDB reformats the data,

performs a calibration of the IR data, adds grid and

orbit/attitude information, and retransmits the data back

scan direction

first pixel

r -

L _

second

pixel

pixel

third pixel

third pixel

Figure 3-1 Oversampling of GOES infrared data

33

to GOES for relay to various user stations. The

retransmitted data contains orbit and attitude parameters

which are required for computing the line and element of

the VISSR data corresponding to a given point on the

surface of the earth.

CHAPTER 4

STUDY AREA AND DATA SETS

Five storm events which occurred over the Kissimmiee

River Basin of Central Florida during the summers of 1979

and 1982 were studied. Figure 4-1 shows the study area and

locations of 32 raingages within this area. Two types of

data were used in this study: (1) hourly rainfall intensity

measured by raingages, and (2) GOES satellite infrared

data.

Hourly Rainfall Data

Hourly rainfall data of the five historic storm events

(Table 4-1) are originally from South Florida Water

Management District (SFWMD). The hourly rainfall data are

accumulated rainfall amounts which occurred between the

nominal hour and the next hour. For example, hourly

rainfall of 9:00 A.M. is the accumulated rainfall that

occurred between 9:00 A.M. and 10:00 A.M.. The geographic

coordinates of raingages inside the study area are listed

in Table 4-2.

GOES Infrared Data

GOES infrared data include satellite images of the

o: Raingage of South Florida Water Management District

(Pixels within the framed area form the training data

set of the pattern classification technique.)

Figure 4-1. The study area and raingage locations.

36

Table 4-1. Historic storm events.

*

DATE DURATION (EST)

8/03/1979 1800 2100

9/09/1982 1400 1800

9/21/1982 1100 1800

9/22/1982 0500 1200

9/26/1982 0500 1000

EST represents Eastern Standard Time.

EVENT INDEX

ST79215

ST82252

ST82264

ST82265

ST82269

Table 4-2. Geographic coordinates

study area.

of raingages inside the

Raingage Index

MRF162

MRF17

MRF14

MRF140

MRF23

MRF168

MRF28

MRF191

MRF155

MRF157

MRF186

MRF185

MRF156

MRF19

MRF161

MRF176

MRF207

MRF35

MRF190

MRF189

MRF143

MRF158

MRF187

MRF159

MRF210

MRF41

MRF160

MRF148

MRF40

MRF188

MRF48

MRF43

Latitude

a I iI

Longitude

0o 1

38

five storm events. IR counts corresponding to brightness

temperatures of pixels were recorded on computer compatible

tapes. Original GOES infrared images have the nominal

resolution of 4*2 nm. To reduce the percentage of overlay,

every other pixel along each scan was resampled to form new

images of 4*4 nm resolution. The Mb enhancement curve

(Figure 2-3) was adapted to show the resampled GOES IR

images in color. Life histories of the five storm events

are shown in Figures 4-2 to 4-6. A computer program was

developed to display digital GOES IR data in picture format

based on the Mb enhancement curve.

Brightness temperatures of pixels can be converted

from IR counts using the following equation:

T = -0.5DN + 56.8 (DN < 176)

T = 144.8 DN (DN > 176 ) (4-1)

where T = brightness temperature (0C),

DN = digital number of IR count.

Pixel Registration of GOES Satellite Data

Each GOES data file consists of IR counts of pixels

along with navigation data of the satellite image itself.

Navigation data is needed for pixel registration. Pixel

registration is the work that associates earth-based

coordinates, usually latitude and longitude, with the

pixels of the image. Navigation data of a satellite image

includes orbital parameters, attitude parameters, spin rate

of satellite, frame geometry, and camera geometry. A

39

computer program developed by SSEC of the University of

Wisconsin was used in this study to convert geographic

coordinates to image coordinates.

40

(a)

Figure 4-2. Life history of Storm-79215. (a) Original

image (b) Enlarged study area.

Figure 4-3. Life history of Storm-82252. (a) Original

image (b) Enlarged study area.

42

(a)

UD)

Figure 4-4. Life history of Storm-82264. (a) Original

image (b) Enlarged study area.

Figure 4-5. Life history of Storm-82265. (a) Original

image (b) Enlarged study area.

(a)

(b)

Figure 4-6. Life history of Storm-82269. (a) Original

image (b) Enlarged study area.

CHAPTER 5

THEORY OF REGIONALIZED VARIABLES--KRIGING/CO-KRIGING

Kriging is a method for optimizing the estimation of a

magnitude which is distributed in space and is measured at

a network of points. Variables that characterize a spatial

phenomenon usually have a behavior much too complex to be

described by usual analytical expressions. Their

variations in space are erratic and often unpredictable

from one point to another. However, their behavior is not

completely random; values taken at neighboring points are

related by a complex set of correlations reflecting the

structure of the underlying phenomenon. Matheron (1971,

1973) first proposed the term 'regionalized variables' to

describe the particular features of these variables.

Let X1, X2,...,Xn be the locations of n points of

measurement. Xi denotes the position vector of point i.

Also let zi=z(Xi) be the value of the variable of interest

Z measured at point i. The problem of point estimation

lies in determining the value of z0 for any point X0 where

measurement is not available. The measured zi values are

considered a realization of a random function Z(X,O). In

order to estimate z0, the structure of the random function

Z(X,0), i.e., its autocovariance function must be

constructed from its realization (zl, z2,..., zn}.

45

46

Second-Order Stationarity and Intrinsic Stationarity

In order to simplify the notation, the realization

variable 0 will be dropped from now on.

If a random function Z(X) is stationary up to order 2

or weakly stationary, then:

1. E[Z(X)]= m, a constant independent of X, and

2. E{[Z(Xi)-m][Z(Xj)-m])= C(Xi-Xj), function of Xi-Xj

only.

For many regionalized variables such as piezometric

head, rainfall intensity, and transmissivity of an aquifer,

the assumption of second-order stationarity is too

restrictive (Delhomme, 1978; Bastin, Lorent, Duque, and

Gevers, 1984; Marsily, 1986). Thus the intrinsic

stationarity is applied to these regionalized variables.

The intrinsic stationarity requires weak stationarity

not for Z(X) itself, but only for its increments. Thus we

need only the hypothesis that, for any vector h,

E[Z(X+h)-Z(X)]= m(h), (5-1)

Var[Z(X+h)-Z(X)]= 27(h) (5-2)

The function j(h) is called the variogram of Z(X).

Assume m(h)=0, i.e. E[Z(X+h)]=E[Z(X)], then

E[Z(X+h)-Z(X)]=0 (5-3)

7(h)=Var[Z(X+h)-Z(X)]/2

= E{[Z(X+h)-Z(X)]2}/2 (5-4)

Computing the Experimental Varioqram

In the case where the mean of Z(X) is a constant, q(h)

47

is defined by Eq. 5-4. To estimate the variogram, we

simply use the measurement points Zi and assume ergodicity

on the increments (i.e., that spatial average can be used

to estimate the realization average).

First we define a certain number of distance classes.

Then, taking all possible pairs of points i and j for each

distance class, we calculate:

1. the number of pairs present in the class,

2. the average distance in the class, and

2

3. the average square increment (Zi-Zj) /2 of

each class.

Thus, from the results of the above steps 2 and 3, the

estimated variogram y(h) can be constructed.

Commonly Used Variogram Models and Their Properties

Just as the autocovariance function must be positive

definite, not any function y(h) can be a variogram. It can

be shown that:

Minus 7 must be conditionally positive definite, i.e.,

for all X1, X2,...,Xn e Rm and for all al,...,an e R, n

n

coefficients satisfying .ai = 0, then

i=1

n n

-E E aiajI(Xi-Xj))0 (5-5)

i=l j=1

Some commonly used variograms are illustrated in Table

5-1 and Figure 5-1.

A variogram, as defined by Eq. 5-4, describes some

spatial characteristics of the random function Z(X): the

48

Table 5-1. Commonly used variogram models.

Model in hA

Spheric model

Variogram 7(h)

whA

w[1.5(h/a)-0.5(h

(h:separating distance)

A<2

/a)3] h

h>a

Exponential model w[l-exp(-h/a)]

Gaussian model w{l-exp[-(h/a)2])

49

3.

2 // 0.5 /

S i -

0 1 2 3 4 0

a)h (b)

1.05- -- ----- 1.0. --r--

0.95 I0.95 -

0 a o 0

h

Cc) (d)

(a) Model in hA

(b) Spherical model

(c) Exponential model

(d) Gaussian model

Figure 5-1. Commonly used variogram models.

50

continuity, the range of influence, and the anisotropies.

The Continuity

The continuity is reflected by the rate of growth of

y(h) for small values of h. Typically, for h=0, -(h)=0

regardless of the variogram. However, variograms may

sometimes exhibit a jump at the origin. This jump at the

origin is called the nugget effect. When the nugget effect

is present, two distinct points Xi and Xj, even though they

are very close to each other, can have significantly

different measurements Z(Xi) and Z(Xj). The nugget effect

may be due to two causes: (1) the micro-regionalization on

a scale much smaller than the spacing of the data points;

and (2) measurement errors (Delhomme, 1978).

The Range of Influence

In some cases, the variogram reaches a limiting value,

called the sill, which is equal to the variance of the data

in the field if the random variable of interest, Z(X), is

second order stationary. The distance at which the sill is

reached is called the range of influence. At distances

beyond the range of influence, correlations between points

are zero.

The Anisotropies

Variograms in different directions may not always be

the same. If they are identical in different directions,

then the continuity of the random function Z(X) is the same

in all directions. Anisotropy can be easily detected by

computing the variogram in different directions. The

51

anisotropy can be eliminated by an appropriate linear

transformation in the coordinate system (David, 1977 and

Clark, 1979).

Parameters Inference of Selected Varioqram Model

Using the experimental variogram in the Kriging method

of estimation may lead to an ill-conditioned coefficient

matrix and, in many cases, completely absurd estimated

values (Bastin and Gevers, 1985). Therefore, it is

necessary to use analytic variogram models, and to infer

the parameters of these models from the data. Bastin and

Gevers (1985) presented several methods for estimation of

variograms. In our study the exponential model of the

following form was selected:

q(h) = w[ 1 exp(-h/a ) ] (5-6)

The exponential model was chosen in this study because it

resembles experimental variograms and its parameter

inferrence is relatively easy. The NLIN procedure (PROC

NLIN) of the Statistical Analysis System (SAS) was used to

infer the least-squares estimates of the exponential

variogram model (SAS, 1985).

Effect of Oversampling of GOES IR Data

on Variogram Estimation

As described in chapter 3, original GOES IR data are

oversampled along the scan. For each satellite image,

every other pixel along each scan is resampled to create a

new image for use in this study. The resampled images have

52

the resolution of 6.87 x 6 km, i.e., the center-to-center

distance of the two adjacent pixels along the scan

direction is 6 km while that distance in N-S direction is

6.87 km. There is 13% overlapped area along the scan

direction. This effect must be considered while

calculating variograms based on image coordinates by

multiplying 0.87 to pixel distances calculated in W-E

direction.

Optimum Estimation by Kriging

We want to find an unbiased estimate Z0* of the

unknown quantity ZO by

n

Z* = EA0iZi (5-7)

i=l

where n is the number of measured points and Aoi is the

weight of Zi when estimating Z0.

Since ZO* is an unbiased estimate, we have

n

EA0i =1 (5-8)

i=l

Also, variance of the estimation error ZO*-ZO can be repre-

sented by the following equation.

Var(Zo* Z0) = E[(ZO* Z0)2]

n n n

= 2Eoyi7(Xi-Xo) oi 0j7(Xi-Xj)

=i=l j=l 01i

(5-9)

Minimization of Var(Zo*-Z0), subject to Eq. 5-8, can

be done by using the Lagrange multipliers. We simply

minimize the expression:

n

E[(Zo*-Z0)2]/2 p(.E AOi-l)

1=1

where p is a new unknown, called a Lagrange multiplier.

The variance of estimation error, subject to the constraint

Eq. 5-8, is minimized when the following equation is

satisfied.

n

E.Aj'Y(Xi-Xj) + P = y(Xo-Xi) i=1, 2,...,n

n

.E Aoi= 1 .

1 A(5-10)

Or, in matrix form:

F1 x A = 2 (5-11)

where

0 '12 113'" in 1

121 0 123 '2n 1

rI =.

7y ... 0 1

nl n2 7 n3 0 1

1 1 1 ... 1 0

A J

01 01

AO I1-

02 02

A= F2 = I

A r2

On On

L '' 1 -

Yij '= (Xi-Xj).

The variance of the error of optimal estimation by

Kriging is:

n

Var(Zo*-ZO) =EAOi(XO-Xi) + u (5-12)

i=1

Once the variogram is known, the unknowns Aqi and p in

Eq. 5-10 or Eq. 5-11 can be solved and Eq. 5-7 can be used

to estimate ZO. Also, by assuming that the distribution of

the estimation error is normal, the confidence interval of

Z0 at a certain level, for example 95% confidence interval

of ZO, can be determined.

Optimum Estimation by Co-Kriqing

Co-Kriging is an estimation technique useful when two

or more correlated variables are measured in the same

field.

Let ZI(X) and Z2(X) be two regionalized variables that

are correlated. The estimation of Z1 by Co-Kriging is

given by a best linear unbiased estimate in the form

n m

ZI*(XO) = ~1AjZ1(Xj) + EA2kZ2(Xk) (5-13)

j=1 k=l

Here, ZI*(Xo) is the estimate of Z1 at location X0.

Z1(Xj) and Z2(Xk) are measured values of Z1 and Z2,

respectively, and A's are Co-Kriging weights.

We now generalize the above equation for N variables

by writing

N n.

Zi*(Xo) = E E AajZa(Xj) (5-14)

a=l j=l

Assuming that the intrinsic stationarity holds for all

Za(X), the Co-Kriging system for estimating Zi(Xo) is given

N n.

E E Akfl (Xj-Xk) + +

a=1 k=l

= 7i(Xo-Xj)

for j=l,...,n and P=1,...,N

(5-15)

and

na 1 if a=i

SAk = if a=i 1,2,...,N (5-16)

k=l 0 if 14i

where yij(h) is called cross-variogram and defined as

Tij(h) = Cov{[Zi(X+h)-Zi(X)],[Zj(X+h)-Zj(X)]}/2

= E{[Zi(X+h)-Zi(X)(X)[Z(X+h)-Zj(X)]/2 (5-17)

(if E[Zi(X)] and E[Zj(X) are both constants.)

The variance of the estimation error is:

Var[Zi*(Xo)-Zi(xo)] =

N n.

E E Auk'ci(XO-Xk) + uPi

a=1 k=1

(5-18)

Myers (1982) proposed the following matrix formula-

tion of Co-Kriging:

n

E 2(Xi-Xj)j + M=2(Xi-Xo)

j=1

i=1,2,...,n

n

Erj = I

j=1

where

(5-19)

n = the number of points where measurements are

available for every variable,

N = the total number of variables,

I = identity matrix of dimension nxn,

0 = zero matrix of dimension nxn,

ri = [Ac]

i

A9a = weight for the ith measured value of

variable p when estimating value of variable

a at X0 (a,=1,2, ...,N).

I(Xi-Xj) = [,EY(xi-Xj)]

A = [Aij], i=l,...,N; j=1,...,N

Or,

2(x-xl -- 2(xljxn) 1

a(Xn-Xi) ... 2(XnfXn) I

I I 0

a(Xi-Xo)

2(Xn-XO)

I

(5-20)

The estimation variance may be written in the

following form:

n

a2 = Tr[ E (XO-Xi)ri ] + Tra

i=l

(5-21)

CHAPTER 6

RAINFALL ESTIMATION--CO-KRIGING APPROACH

This chapter describes detailed procedures of the

proposed scheme for rainfall estimation using GOES

satellite infrared data and raingage measurements. The

estimation scheme consists of two phases. Phase I

partitions GOES pixels into two classes, precipitation-free

and precipitating pixels, using an infrared temperature

threshold and a pattern classification technique. Phase II

estimates pixel-averaged rainfall (PXLRN) using the Co-

Kriging method. Two random fields, PXLRN and the coldest

cloud-top temperature (CCTT), are involved in Co-Kriging.

Several criteria used to verify the model validity and to

evaluate estimation results are also described in this

chapter.

Phase I--Pixels Classification

Delineating Cloud-covered Area

There are three major features in each GOES satellite

image: clouds, land surface, and sea surface. Pixels

representing land surface and sea surface are apparently

precipitation-free. Since cloud pixels are generally

colder than pixels of the other two features, GOES IR

temperature can be used to delineate cloud-covered areas.

57

58

As described in chapter 2, threshold temperature of -200C

was chosen by Griffith et al. (1981) and Weiss and Smith

(1987). Threshold temperatures of -320C and -15.50C were

also used by Scofield and Oliver (1977) and Wylie (1979).

These threshold values were chosen because the cloud areas

measured from them exhibited time-variations that closely

correlated with the radar measured rain volumes. Wylie

(1979) indicated that the cloud areas and consequently rain

rates are extremely sensitive to the colder thresholds

while at warmer thresholds the sensitivity was reduced. In

this study -200C was adopted as the threshold. Pixels with

IR temperatures higher than the threshold value are

considered precipitation-free. However, not all pixels

with IR temperatures lower than -200C are precipitable. A

statistical pattern classification technique is needed in

order to further classify these pixels into two categories:

precipitation-free and precipitating.

Determining Classification Features

Features used in the pattern classification technique

were derived from pixel groups as shown in Figure 6-1.

According to the results of Wu et al., (1985) the coldest

cloud-top temperature within a pixel group best separates

precipitation-free and precipitating pixels. Similar

results were found by Weiss and Smith (1987). In addition

to CCTT, both methods require other features in the pattern

classification technique. IR features [I], [4], [10],

First

pixel group

< ------

central

pixel

Second

pixel group

a pixel group

GOES satellite image

Figure 6-1. A pixel group used to determine values of

classification features.

I I I I I

60

[14], and [15] (Table 2-2) were also used by Wu et al.

(1985). Weiss and Smith included the average cloud-top

temperature and its change over a period of 30 minutes in

their classification scheme. However, texture features of

Wu et al. (1985) won't be reliable if GOES IR data with

6.87*6.87 km resolution are to be used. Within each 20*20

km grid there are only 9 pixels and it is inadequate to

derive texture features based on the histogram of radiance

differences of the pixel group. Also, the 0.4 inch/hr

cutoff rainrate of Weiss-Smith method is not practical from

the hydrologists point of view.

One major difference between raingage measurements and

inferred rainfall using satellite data is their represen-

tative areas. Raingage measurements are actually point

measurements while satellite inferred rainfalls are

averages over an area equivalent to the IFOV of satellite

data. In order to verify the satellite-based rainfall

estimation scheme, average rainfall over a pixel size must

be obtained using raingage measurements. PXLRN calculated

from raingage measurements are considered ground-truths in

this study since their true values are never known.

Consider a group of 9 pixels as shown in Figure 6-1.

PXLRN of the central pixel, pixel 0, is determined as the

arithmetic average of all raingage measurements within the

pixel group. Two other features, CCTT and standard

deviation (STD_DEV) of cloud-top temperatures of pixels

within a pixel group, are calculated for all gaged pixels,

61

i.e., pixels including at least one raingage in their

coverages.

In this study, features CCTT and STD_DEV are used in

the statistical pattern classification technique to

separate precipitation-free and precipitating pixels. CCTT

detects the region of active convection where rainfall is

most likely to occur (Griffith et al., 1978 and Weiss and

Smith, 1987). STD_DEV developed in this study is a measure

of the local variation of IR cloud-top temperature in a

pixel group which can be used to detect the

expanding/contracting status of cloud systems.

CCTT and STD_DEV are calculated for all gaged pixels

of images of the five storm events studied. However, only

pixels inside the framed area, the training area, as shown

in Figure 4-1 are used to form the training data set of the

pattern classification technique. The rest of gaged pixels

are used to test the classifier. Since the raingage

density of the training area is much higher than that of

the test area, its PXLRN are considered more reliable.

Mathematical Concepts of Pattern Classification

Pattern classification is also called statistical

discriminant analysis. It is essentially a decision-making

process with data that can exhibit considerable statistical

variability. The decision made at each satellite pixel (in

our case) must minimize some kind of error criterion

throughout the classified area, i.e., over a large number

62

of individual pixel classifications.

An intuitively satisfying and mathematically tractable

classification theory having the above property is maximum

likelihood, or Bayes optimal, classification. The

mathematics of Bayesian classification is briefly described

below.

Suppose n features of a GOES satellite image are

measured and we want to decide to which of the two classes

(for example, precipitation-free and precipitating) a pixel

belongs.

Let p(Xli) represent the conditional probability

density function of feature vector X of class i and p(i)

represent the probability that class i occurs in the image

area of interest. p(i) is also called the a priori

probability of class i. To make a classification decision

for a pixel, we need to know the a posteriori probabilities

that the pixel belongs to each of the two classes, given

that the pixel has the feature value X.

According to the Bayes rule, the a posteriori

probability is:

p(iIX) = p(Xli)p(i)/p(X) (6-1)

where

2

p(X) = E p(Xli)p(i) (6-2)

i=l

Since p(X) is the same for both classes in Eq. 6-1, we

then define a discriminant function Di(X):

Di(X) = ln[p(Xli)p(i)] (6-3)

Based on the discriminant function Di(X), the Bayes

63

decision rule can be stated as

a pixel belongs to class 1 if D1(X) > D2(X)

a pixel belongs to class 2 if D1(X) < D2(X).

It can be shown (Duda and Hart, 1973) that the Bayes

decision rule minimizes the average probability of error

over the entire classified data set, if all the classes

have normal probability density functions.

In the case that p(Xli) is normal; i.e.,

1

p(X|i) = ----------------exp[-1/2(X-Mi)t Ei-(X-Mi)]

Ei 11/2(2r) k/2

(6-4)

where

X = pixel feature vector,

Mi = mean vector of X for class i,

Ei = covariance matrix of X for class i and

k = number of features,

the Bayes optimal discriminant function becomes

(Schowengerd, 1983)

k 1

Di(X) = n[p(i)] -ln(27r) -IlnIEi

2 2

(X-Mi) Ei-1(X-Mi)

2

(6-5)

Classification of Cloud Pixels Using Bayes Optimal

Classifier

Prior to the performance of the pattern classification

technique, several cutoff rainrates were chosen to define

64

precipitation-free and precipitating pixels. Rainrates of

0.00, 0.01, 0.02, 0.03, and 0.04 inch/hr were tested.

Those pixels with PXLRN less than the chosen cutoff

rainrate were considered precipitation-free. Bayes optimal

classifier was used to perform classifications and a priori

probabilities of precipitation-free and precipitating

pixels were both set equal to percentages of their

occurrences in all images used for classification.

The classification results for both training and test

data sets are shown in Table 6-1. Misclassification rates

of precipitation-free pixels vary from 71% to 100% for the

training data set and from 80% to 100% for the test data

set. Misclassification rates for precipitating pixels vary

from 19% to 33% for the training data set and from 17% to

35% for the test data set. Overall misclassification rates

vary from 34% to 52% for the training data set and from

34% to 56% for the test data set.

Average values of CCTT, and STD_DEV over the whole

study area were calculated for each satellite image in

order to study the time evolution of these parameters and

to evaluate the effects they may have on pattern

classification. Average STD_DEV reveals a very significant

pattern (Figure 6-2). It can be seen clearly that average

STD_DEV decreases as the cloud system diminishes. Though

average STD_DEV of Storm-82269 increased from 0500EST to

0800EST and then decreased, their values are relatively

small. Actually, Storm-82269 formed before 0000EST and

65

Table 6-1. Misclassification rates of pattern classification

technique-all images

Classification features: CCIT and STD DEV

1

Training Data Set (N=361)

Cutoff rainrate

(inch/hr)

0.00

0.01

0.02

0.03

0.04

Classes which pixels come from

Precipitation-

free Precipitating

3 4

0.89 (0.21) 0.19 (0.79)

0.91 (0.28) 0.22 (0.72)

0.92 (0.35) 0.24 (0.65)

1.00 (0.41) 0.19 (0.59)

0.71 (0.47) 0.33 (0.53)

Test Data Set (N=395)

Cutoff rainrate

(inch/hr)

0.00

0.01

Classes which pix(

Precipitation-

free

0.96

0.97

els come from

Precipitating Overall

0.18 0.34

0.19 0.41

0.02 0.93 0.22 0.47

0.03 1.00 0.17 0.51

0.04 0.80 0.35 0.56

1Number of pixels used for classification.

2Misclassification rates weighted by a priori probabilities.

3A priori probability of precipitation-free pixels.

4A priori probability of precipitating pixels.

2

Overall

0.34

0.41

0.48

0.52

0.51

~-- -I-------------I

c3B--e-Os

ST79215

ST82252

ST82264

ST82265

ST82269

4.00 6.00 8.00 10.00 12.00 14.00 16.00

AVERAGE STD-DEV (C)

Figure 6-2. Evolution of average STD_DEV.

2400

2000-

1600:

1200

800-

400

7'

I I

I ~

0.00

2.00

I

'''''l'''''''

67

started diminishing before 0500EST. Unfortunately, GOES

satellite images before 0500EST were not available.

The evolution of average STD_DEV suggests that the

life history of a cloud system should be taken into

consideration to separate cloud pixels into different

classes. By considering the values of average STD_DEV and

examining the cloud coverages of the study area, we then

divided the total of 31 GOES satellite images into two

groups: expanding images and contracting images. Expanding

clouds have higher IR temperature variation than do

contracting clouds. Table 6-2 lists the expanding/con-

tracting status of each GOES satellite image.

For the training data set, histograms of PXLRN of

expanding and contracting images (Figure 6-3) indicate that

approximately 90% of PXLRN of contracting images are less

than 0.2 inch/hr; only 60% of PXLRN of expanding images are

less than 0.2 inch/hr. Most heavy rainfalls as defined by

0.2 inch/hr occurred during expanding stages. Some

statistics of the two histograms are listed in Table 6-3.

The t-test for the equality of the two means rejects the

null hypothesis at the significant level a=0.05, i.e.,

PXLRN averages of expanding and contracting images are

significantly different.

Pattern classification technique along with different

cutoff rainrates was then applied to images of expanding

and contracting stages respectively. Comparing to the

results of previous classification (Table 6-1),

Table 6-2.

Expanding/contracting status of GOES satellite

images

Storm event

79215

Time of recording

1800EST

1900EST

2000EST

2100EST

Status

expanding

expanding

contracting

contracting

82252 1400EST expanding

1500EST contracting

1600EST contracting

1700EST contracting

82264 1100EST expanding

1200EST expanding

1300EST expanding

1400EST expanding

1500EST contracting

1600EST contracting

1700EST contracting

1800EST contracting

82265 0500EST expanding

0600EST expanding

0700EST contracting

0800EST contracting

0900EST contracting

1000EST contracting

1100EST contracting

1200EST contracting

0500EST

0600EST

0700EST

0800EST

0900EST

1000EST

contracting

contracting

contracting

contracting

contracting

contracting

82269

100

0

z

CONTRACTING

EXPANDING

S 60

W o

4

CE

40

H

2 20

:3

l-

20

II l 1I1 I I I I II i l 'I

0 0.2 0.4 0.6 0.6

PXLRN (INCH/HR)

Figure 6-3. Histograms of PXLRN of expanding and contracting images.

1. 2

70

Table 6-3. Statistics of PXIRN of expanding and contracting

images

Expanding Contracting

Number of pixels 81 280

Average (inch/hr) 0.22 0.08

Variance (inch2/hr2) 0.05 0.01

Median (inch/hr) 0.15 0.05

71

misclassification rates were significantly reduced (Tables

6-4 and 6-5). For expanding images, misclassification

rates of precipitation-free pixels vary from 35% to 76% for

the training data set and from 65% to 86% for the test data

set. Misclassification rates of precipitating pixels vary

from 8% to 23% for the training data set and from 13% to

30% for the test data set. Overall misclassification rates

vary from 27% to 32% and from 34% to 44% for training data

set and test data set, respectively. For contracting

images, misclassification rates of precipitation-free

pixels vary from 39% to 89% for the training data set and

from 63% to 91% for the test data set. Misclassification

rates of precipitating pixels vary from 2% to 24% for the

training data set and from 4% to 26% for the test data set.

Overall misclassification rates vary from 18% to 32% and

from 21% to 44% for training and test data sets,

respectively. In general, misclassification rates are

higher for test data set than for training data set. This

is not unexpected considering that raingage density of the

test area is lower than that of the training area. If more

raingage measurements were taken in the test area, some

zero-PXLRN pixels may in fact have non-zero PXLRN.

Both image categories consist of many zero-PXLRN

pixels. These pixels represent area of non-occurrence

while the rest represent area of occurrence. In order to

study the spatial variation of rainfall, zero-PXLRN pixels

must be identified. More than 70% of pixels in the

Table 6-4. Misclassification rates of pattern classification

technique-expanding images

Classification features: CCIT and STD DEV

Training Data Set (N=81)

Cutoff rainrate

(inch/hr)

0.00

0.01

0.02

0.03

0.04

Cutoff rainrate

(inch/hr)

0.00

0.01

0.02

0.03

0.04

Classes which pixels ccme from

Precipitation-

free Precipitating

0.76 (0.29) 0.08 (0.71)

0.76 (0.31) 0.13 (0.69)

0.41 (0.35) 0.19 (0.65)

0.35 (0.38) 0.23 (0.62)

0.38 (0.39) 0.20 (0.61)

Test Data Set (N=95)

Classes which pixels come from

Precipitation-

free Precipitating

0.85 0.13

0.86 0.19

0.80 0.20

0.65 0.29

0.67 0.30

Overall

0.27

0.32

0.27

0.28

0.27

Overall

0.34

0.40

0.41

0.43

0.44

- --- --"--

Table 6-5. Misclassification rates of pattern classification

technique-contracting images

Classification features: CCIT and STD DEV

Training Data Set (N=280)

Cutoff rainrate

(inchhr)

0.00

0.01

0.02

0.03

0.04

Cutoff rainrate

(inch/hr)

0.00

0.01

0.02

0.03

0.04

Classes which pixels come from

Precipitation-

free

0.89

0.71

0.64

0.57

0.39

(0.19)

(0.27)

(0.36)

(0.42)

(0.49)

Precipitating

(0.81)

(0.73)

(0.64)

(0.58)

(0.51)

0.02

0.03

0.09

0.15

0.24

Test Data Set (N=300)

Classes which pixels come from

Precipitation-

free Precipitating

0.91 0.04

0.88 0.10

0.81 0.11

0.75 0.17

0.63 0.26

Overall

0.18

0.21

0.29

0.32

0.32

Overall

0.21

0.31

0.36

0.41

0.44

--

74

training data set were correctly classified for both image

categories if the cutoff rainrate of 0.02 inch/hr was used.

The 0.02 inch/hr cutoff rainrate is in accordance with the

value of operational radar precipitation levels: no rain (0

< R < 0.02 inch/hr); light rain (0.02 < R < 0.2 inch/hr);

and heavy rain (0.2 inch/hr < R) (Wu et al., 1985). Thus,

cutoff rainrate of 0.02 inch/hr is used to identify

precipitation-free and precipitating pixels.

A conceptual model of rainfall variation over space is

proposed as follows:

The PXLRN rainfall field consists of two subsets, zero

subset and non-zero subset, based on the cutoff rainrate of

0.02 inch/hr. The two subsets are characterized by

classification features CCTT and STDDEV. Zero subset

represents area of non-occurrence and non-zero subset

represents area of occurrence. Spatial variation of PXLRN

of non-zero subset is independent of PXLRN of zero subset

and zero subset plays no role in the estimation of non-zero

PXLRN.

Thus, estimation of PXLRN of ungaged pixels can be

made by the following steps.

1. Calculate values of the two features, CCTT and

STD_DEV, for pixels of interest.

2. Perform pattern classification on these pixels

to identify zero PXLRN pixels and non-zero PXLRN

pixels.

3. Assign zero PXLRN to identified zero PXLRN pixels.

75

4. Estimate PXLRN for identified non-zero PXLRN pixels

using the Co-Kriging method described in Phase II.

Figure 6-4 illustrates that cloud-covered area is

delineated using IR temperature threshold and then

precipitating area is defined by features CCTT and STD_DEV.

Phase II--Estimation of PXLRN for Ungaged Pixels

Using Co-Kriging Method

Formulation of Co-Kriging Equation for PXLRN Estimation

As described in Phase I, CCTT and STD_DEV are used in

the pattern classification technique to delineate the

precipitating area. Within the precipitating area, PXLRN,

CCTT, and STD_DEV are available for gaged pixels. However,

for ungaged pixels only the last two features are available

and PXLRN must be estimated.

Assume PXLRN and CCTT are two random variables and

their variations over space form two random fields: PXLRN

field and CCTT field. Our goal is to estimate PXLRN of

ungaged pixels based on the spatial autocorrelation of

PXLRN and CCTT and spatial cross-correlation of the two.

The following Co-Kriging equation is employed for PXLRN

estimation:

n

R0(t) = E (A R(t)R.(t) + A ,(t)T i(t))

1=1 Ri i Ti i

+ T (t)T (t)

TO 0

(6-6)

ZERO PXLRN REGION

NON-ZERO REGION

Boundary of precipitating

area defined by OCTT and

STD DEV.

Boundary of cloud-covered

-area defined by IR

temperature threshold.

Figure 6-4, A cloud system shows precipitation-free

and precipitating areas.

where

R0(t) = PXLRN of ungaged pixel at location X0;

Ri(t) = PXLRN of gaged pixel at location Xi;

Ti(t) = CCTT of gaged pixel at location Xi;

T0(t) = CCTT of ungaged pixel at location X0;

ARi(t) = Co-Kriging weight of Ri(t);

ATi(t) = Co-Kriging weight of Ti(t);

n = number of gaged pixels involved;

X = position vector; and

t = time index.

The above equation can also be represented as:

n

RO (t) = [A (t)R.(t) + A (t)T .(t)] (6-7)

i=0 Ri 1 Ti 1

with A R(t) =0.

Since Co-Kriging estimate is optimal and unbiased,

variance of estimation error, Var[RO(t)-RO(t)], must be

minimized subject to the following constraints:

n

EA Ri(t) = 1, and (6-8)

i=0 Rl

n

E A (t) = 0 (6-9)

i=0 Ti

Under the assumptions that R(t) and T(t) are both

intrinsic stationary random fields and have constant means,

the Co-Kriging weights ARi and ATi can be obtained by

solving the Co-Kriging equation system:

n n

E A Ri RR(Xk-Xi) + E A Ti (Xk-Xi)

i=0 Ri RR i=0 Ti RT

+ p R RR(XO-Xk)

n n

E A (Xk-Xi) + E A ,y T (Xk-Xi)

i=0 Ri TR i=0 Ti TT

+ T = IRT(XO-Xk)

(k=0,l,...,n)

n

E A = 1

i=0 Ri

n

E A = 0

i=0 Ti

(6-10)

where

RR = variogram of PXLRN field;

TT = variogram of CCTT field;

TT

T = cross-variogram of PXLRN and CCTT ; and

RT

(7RT TR )

R T= Lagrange multipliers.

The time variable t does not appear in the above

equation system since it must be solved for each individual

image.

Random Field Analysis

Before Eq. 6-10 can be solved, variograms RR and

TT and cross-variogram RTmust be obtained. Variograms

TT RT

79

characterize the spatial variabilities of random fields.

Non-zero PXLRN and CCTT of gaged pixels are used to

construct variograms and cross-variogram.

Estimation of variograms. By definition, q(h) =

1

- Var[Z(X+h) Z(X)]. Under the assumption that E[Z(X)] is

2 1

a constant for any X, 7(h) = E[Z(X+h) Z(X)]2 For

2

each satellite image, variograms of PXLRN and CCTT can be

obtained by following steps described in chapter 5. Figure

6-5 shows the experimental variograms of PXLRN of Storm-

82269. Experimental variograms are not, in general,

positive definite functions (Armstrong and Jabin, 1981).

Certain types of theoretical variogram models as described

in chapter 5 must be adopted and experimental variograms be

fitted.

Here, a problem must be addressed before variogram

fitting is performed. Of all gaged pixels, some of them

may have zero PXLRN and are not to be included for

variogram fitting. For images including only a few non-

zero PXLRN gaged pixels, (for example less than 10 pixels)

experimental variograms thus derived can not be considered

reliable.

Variograms, as shown in Figure 6-5, are functions of

separating distance and time as well. For those satellite

images which reliable variograms can not be obtained,

adopting variograms of previous images is not an adequate

alternative because variograms are not stable over time.

Dimensionless variograms. The variogram of Eq. 5-6

0.03

. 0.02

0.01

0.00

h (pixel)

Figure 6-5. Experimental variograms of PXLRN of Storm-82269.

81

has two parameters, w and a. During the life history of a

cloud system, variograms vary from time to time and can be

represented as

7(h,t) = w(t)[l-exp(-h/a(t))] (6-11)

where t is the time variable. Parameters w(t) and a(t) are

both functions of t. Changes of variograms are due to

changes of w(t) and a(t) over time. Figures 6-6 and 6-7

illustrate changes of an exponential variogram when w and a

are independent of t, respectively. Figure 6-8 shows the

changes when both w and a are functions of t.

Changes of a(t) do not change the sill of a variogram.

It only changes the speed at which the variogram approaches

its sill. On the contrary, changes of w(t) do not change

the approaching speed, but change the sill of a variogram.

It can be seen clearly in Figure 6-8 that changes of

variograms due to changes of a(t) are relatively small

compared to that due to changes of w(t).

Standard deviation and average of PXLRN of each GOES

image were calculated. Average PXLRN and standard

deviation of PXLRN are highly correlated. Figure 6-9

indicates that large standard deviation, and hence

variance, of PXLRN comes with large average PXLRN. This is

a clear indication that a variogram is affected by the

spatial autocorrelation as well as the variance of the

random field. Therefore, we attribute all time-dependent

variations of a variogram to w(t) and Eq. 6-11 becomes

7(h,t) = w(t)[l-exp(-h/a)] (6-12)

1.20

0.80

-c

L-

0.40

0.00 -

0.00

Figure 6-6.

Changes of the

independent of

12.00

Distance h

exponential variogram when w is

time.

4.00

3.00

S2.00

1.00

0.00

Figure 6-7.

Distance h

Changes of the exponential variogram when a is

independent of time.

2.50

2.00

1.50

1.00

0.50

0.00

Figure 6-8.

) 4.00 6.00 8.00 10.00

Distance h

Changes of the exponential variogram when both w and a

are functions of time.

w=2

w = 1

w = 0.5

12.00

0.20 0.30

AVERAGE OF PXLRN (INCH/HR)

Figure 6-9. Average of PXLRN versus standard deviation of PXLRN.

0.60

rY

i

I 0.50

z

z

' 0.40

x

0-

a

I-

0

z 0.30

0

o 0.20

z

<0.10

0.00

0.00

86

Parameter w(t) changes from image to image and is the

parameter which takes into consideration the variance of

the random field. Parameter a takes into consideration the

spatial autocorrelation and is a constant during the life

history of a cloud system.

For the case that Z(X,t) is second order stationary,

q(h,t) = C(0,t) C(h,t) (6-13)

Var[Z(X,t)] = az2(X,t) = y(oo,t)

= w(t) (6-14)

where C(h,t) = Cov[Z(X,t),Z(X+h,t)]. If standardized

values of Z(X,t), i.e., Z(X,t)/cz(X,t), are used to

construct the variogram, then the sill of the new variogram

equals unity, no matter which satellite images are used.

In other words, there exists a unique dimensionless

variogram, if Z(X,t) is second order stationary and the

variogram model of Eq. 6-12 is adopted.

PXLRN and CCTT were assumed second order stationary in

this study since the orographic influence is not

significant in Florida's flat area. Using dimensionless

variogram has two advantages. Firstly, dimensionless

experimental variograms of different satellite images can

be considered realizations of a single population. Thus,

PXLRN or CCTT of different images can be combined to

estimate theoretical variograms. This would help solve the

problem of not having enough numbers of non-zero PXLRN

gaged pixels for each individual satellite image and the

estimated variograms would also be more reliable.

87

Secondly, since the dimensionless variogram is unique,

variogram fitting only needs to be performed once and

original variograms can be easily converted from the

dimensionless variogram. The following section describes

how original variograms can be converted.

Let Z(X,t) be a second order stationary random field.

Thus, E[Z(X,t)] and Var[Z(X,t)] are both independent of X.

Original variogram of Z(X,t) can be represented as

1

y(h,t) = E[Z(X+h,t)-Z(X,t)]2

2

= w(t)[l-exp(-h/a)] (6-15)

Define Z*(X,t)-Z(X,t)/oaz(t) where az2(t) is the

variance of Z(X,t). Dimensionless variogram can then be

expressed as

y*(h,t) = 1 E[Z*(X+h,t)-Z*(X,t)]2

1 1

1 ---- { E[Z(X+h,t)-Z(X,t)]2 )

az2(t) 2

1

S------ (h,t) (6-16)

az2 (t)

S *(h,t) =------w(t)[l-exp(-h/a)]

az2(t)

= w (t)[l-exp(-h/a)] (6-17)

( w(t) = (t)z2(t) )

The above derivation shows that original variograms

can be obtained by simply multiplying the variance of

88

Z(X,t), az2(t), to the dimensionless variogram.

Suppose a cloud system lasts for d hours and there are

ni non-zero PXLRN gaged pixels in the ith hour satellite

image. Let Z(X,t) be the random field of PXLRN or CCTT.

We then estimate the mean and variance of Z(X,t) by the

following equations:

^ 1 ni

E[Z(X,i)] = --- Z(Xj,i)

ni j=1

A

= pz(i) (6-18)

A 1 ni A

Var[Z(X,i)] = --- E [Z(Xj,i) Az(i)]2

ni j=l

A

= A2(i) i= ,2,...,d (6-19)

For the ith satellite image there are ni(ni-l)/2 pairs

of non-zero PXLRN gaged pixels [Z(Xp,i), Z(Xq,i)], p#q. A

d

total of ni(ni-l)/2 = N pairs are then used to

construct the dimensionless variogram of Z(X,t) using the

method of nonlinear fitting.

Figures 6-10 to 6-19 demonstrate dimensionless experi-

mental variograms of PXLRN and CCTT of the five storm

events used in this study. PROC NLIN of SAS was used to

infer parameters of dimensionless variograms directly from

pixel pairs. Parameter values of fitted dimensionless

variograms of PXLRN and CCTT of different storm events are

listed in Table 6-6. Figures 6-20 and 6-21 present fitted

dimensionless variograms of PXLRN and CCTT respectively.

Dimensionless variograms of PXLRN show a very consistent

pattern among different storms. This is an indication