30 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 22

Constructing Three-Dimensional Multiple-Radar Reflectivity Mosaics:

Examples of Convective Storms and Stratiform Rain Echoes
JIAN ZHANG
Cooperative Institute for Mesoscale Meteorological Studies, University of Oklahoma, Norman, Oklahoma

KENNETH HOWARD
National Severe Storms Laboratory, Norman, Oklahoma

J. J. GOURLEY
Cooperative Institute for Mesoscale Meteorological Studies, University of Oklahoma, Norman, Oklahoma

(Manuscript received 2 April 2004, in final form 20 July 2004)

ABSTRACT
The advent of Internet-2 and effective data compression techniques facilitates the economic transmission
of base-level radar data from the Weather Surveillance Radar-1988 Doppler (WSR-88D) network to users
in real time. The native radar spherical coordinate system and large volume of data make the radar data
processing a nontrivial task, especially when data from several radars are required to produce composite
radar products. This paper investigates several approaches to remapping and combining multiple-radar
reflectivity fields onto a unified 3D Cartesian grid with high spatial (ⱕ1 km) and temporal (ⱕ5 min)
resolutions. The purpose of the study is to find an analysis approach that retains physical characteristics of
the raw reflectivity data with minimum smoothing or introduction of analysis artifacts. Moreover, the
approach needs to be highly efficient computationally for potential operational applications. The appro-
priate analysis can provide users with high-resolution reflectivity data that preserve the important features
of the raw data, but in a manageable size with the advantage of a Cartesian coordinate system.
Various interpolation schemes were evaluated and the results are presented here. It was found that a
scheme combining a nearest-neighbor mapping on the range and azimuth plane and a linear interpolation
in the elevation direction provides an efficient analysis scheme that retains high-resolution structure com-
parable to the raw data. A vertical interpolation is suited for analyses of convective-type echoes, while
vertical and horizontal interpolations are needed for analyses of stratiform echoes, especially when large
vertical reflectivity gradients exist. An automated brightband identification scheme is used to recognize
stratiform echoes. When mosaicking multiple radars onto a common grid, a distance-weighted mean scheme
can smooth possible discontinuities among radars due to calibration differences and can provide spatially
consistent reflectivity mosaics. These schemes are computationally efficient due to their mathematical
simplicity. Therefore, the 3D multiradar mosaic scheme can serve as a good candidate for providing
high-spatial- and high-temporal-resolution base-level radar data in a Cartesian framework in real time.

1. Introduction of the WSR-88Ds, it becomes inconvenient when

coanalyzing radar data with other observational data
The deployment of the Weather Surveillance Radar- sources (e.g., satellite data) or when combining data
1988 Doppler (WSR-88D; http://www.roc.noaa.gov/) from multiple radars. Numerous approaches and tech-
network has provided meteorologists critical informa- niques have been developed for gridding WSR-88D
tion toward the issuance of warnings for tornadoes, se- and other ground-based or airborne radars’ observa-
vere storms, and flash floods. Spherical coordinates are tions in Cartesian space. The commonly used interpo-
most commonly used to store the raw radar observa- lation techniques include 1) nearest neighbor (e.g.,
tions. While a spherical coordinate system accommo- Jorgensen et al. 1983), 2) linear interpolation (e.g.,
dates the data collection method (i.e., scan strategies) Fulton 1998; Mohr and Vaughn 1979; Miller et al.
1986), 3) the Cressman weighting scheme (e.g., Wey-
Corresponding author address: Jian Zhang, NSSL, 1313 Halley gandt et al. 2002a), and 4) the Barnes or exponential
Circle, Norman, OK 73069. weighting scheme (e.g., Shapiro et al. 2003; Askelson et
E-mail: [email protected] al. 2000). Software packages, such as the Sorted Posi-

© 2005 American Meteorological Society

JTECH1689
JANUARY 2005 ZHANG ET AL. 31

tion Radar Interpolator (SPRINT; Miller et al. 1986) adequately monitor the evolution of storm and associ-
and REORDER (http://www.atd.ucar.edu/rdp/home/ ated characteristics. The Collaborative Radar Acquisi-
reorder.html), that were developed at the National tion Field Test (CRAFT; Droegemeier et al. 2002) proj-
Center for Atmospheric Research included one or ect successfully demonstrated the feasibility of trans-
more of the aforementioned analysis schemes and mitting full-volume scan data from all radars in the
provide many researchers (e.g., Bluestein and Gaddy NWS WSR-88D network to a center facility economi-
2001; Ziegler et al. 2001; Wurman and Gill 2000) tools cally and in real time. This provides the ability to inte-
for radar data analysis. grate the full-resolution base-level data from multiple
Smoothing and filtering characteristics of these radars onto a common 3D framework. The 3D mosaic
simple schemes are studied in Trapp and Doswell grid can benefit forecasters, meteorologists, and re-
(2000) and Askelson et al. (2000) using theoretical ap- searchers with a wide variety of products and displays,
proximations and observation system simulation ex- including flexible horizontal or vertical cross sections in
periments. More sophisticated analysis schemes include addition to regional rainfall maps. High-resolution re-
statistical (Heymsfield 1976) and variational (Gao et al. flectivity analyses can also serve as an important source
1999) approaches. A majority of the Cartesian interpo- in data assimilations for convective-scale numerical
lation schemes were developed for specific research ap- weather modeling (e.g., Zhang 1999; Weygandt et al.
plications whereby the parameters in the analysis 2002b) over large domains and for merging conven-
schemes are usually dependent on the applications’ ob- tional datasets (e.g., lightning strike information, objec-
jective. For example, some degree of smoothing in tively analyzed rawinsonde observations, numerical
analyses is necessary for reducing errors in single- forecast simulation fields, etc.) with the radar data.
A radar reflectivity observation is not simply a point
Doppler wind retrievals (Zhang and Gal-Chen 1996;
observation. It is an integrated electromagnetic power
Gao et al. 1999). This type of smoothing, however, may
return from scatters in a radar sample volume (Rine-
be undesirable for severe storm applications, such as
hart 2001) or resolution volume (Doviak and Zrnic
storm cell tracking. The purpose of this study is to find
1993). The sample volume is also referred to as a radar
a scheme that can generate a radar reflectivity analysis
data bin in this paper. Due to the spherical geometry of
on a three-dimensional (3D) high-resolution Cartesian
radar sampling, the size of radar sample volumes in-
grid that satisfies the following criteria. creases with increasing range while the sample resolu-
1) The gridded reflectivity data must retain, as much as tion decreases with increasing range. Figure 1 shows the
possible, the important convective-scale storm struc- range height coverage for the WSR-88D volume cov-
tures evident in the raw radar data with minimal erage pattern (VCP) 21 (OFCM 2003) for standard at-
smoothing. mospheric refractive conditions. Within ⬃20 km of the
2) The analyses should be physically realistic with radar, centers of radar data bins are less than 1 km
minimal induced artifacts. apart, while at ranges of ⬃150 km, data are spaced
3) The scheme must be computationally efficient and more than 50 km horizontally (between points A and
economic for real-time operational applications. A⬘) and more than 2.5 km vertically (between B and
B⬘). This nonuniformity in data spacing makes the
With these criteria, the analyses should provide end
users with high-resolution radar reflectivity data fields
that are comparable to the raw data with the advantage
of a Cartesian coordinate system. A Cartesian coordi-
nate system provides a common framework in which
other observational datasets can be merged and cross-
correlated. This facilitates the creation of multisensor
algorithms and applications that use the strength of in-
dividual observations collectively to provide more
physically and scientifically sound depictions of meteo-
rological phenomena than a single observing system.
The transformation of radar data from a spherical
coordinate to a Cartesian grid provides a more direct
approach of combining multiple radars onto a common
grid. Forecasters are often responsible for areas that
encompass multiple-radar umbrellas [e.g., the National
Weather Service (NWS) County Warning Areas and
the Federal Aviation Administration (FAA) Air Route
Traffic Control Centers]. Moreover, the life cycle of an FIG. 1. WSR-88D sampling distributions for VCP 21 assuming a
individual storm or storm system may span a region beam propagation model under standard atmospheric refraction
that requires observations from two or more radars to conditions.
32 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 22

choice of an interpolation scheme and associated filter modes, each of which has prespecified elevation angles
somewhat dependent on the objective of the applica- (Table 1). Figure 2 shows vertical cross sections of
tion. Trapp and Doswell (2000) evaluated the error beam propagation paths and power density distribu-
characteristics of nearest-neighbor, Barnes-type, and tions for the VCP11 and VCP21 scan modes. The
Cressman-type interpolation schemes using simulated propagation of the beams is assumed to follow the 4/3-
radar fields. Their results show that a nearest neighbor effective earth radius model (Doviak and Zrnic 1993):
scheme gives the smallest root-mean-square (rms) er-
rors between the analysis and the observations, but the 4
ae ⫽ a, 共2兲
spatial scales of the error fields are nonuniform. A 3
heavy Barnes filter based on the poorest resolution of
the reflectivity data results in an error field of uniform h ⫽ 共r2 ⫹ a2e ⫹ 2rae sin␪e兲1Ⲑ2 ⫺ ae, 共3兲

冉 冊
spatial scales but with significantly less high-resolution
information (Askelson et al. 2000) than in the raw data. r cos␪e
s ⫽ ae sin⫺1 , 共4兲
This is not desirable for severe storm tracking applica- ae ⫹ h
tions, nor for convective-scale numerical weather as-
similation applications because higher-resolution infor- where a represents the earth’s radius, ae represents the
mation available at shorter ranges is discarded. 4/3-effective earth radius, h is the height of the center of
This paper explores various objective analysis meth- the radar beam, and s is the distance between the radar
ods applied to WSR-88D data across a spectrum of and the projection of the bin along the earth surface.
weather regimes with comparisons between nearest- All the computations in this paper assume this standard
neighbor and two linear interpolation schemes. Results atmospheric beam propagation model.
from several mosaicking techniques, including nearest Equation (1) shows that the radar power peaks along
neighbor, maximum value, and distance-weighted the center of the beam and decreases in the azimuthal
mean, are presented for comparison. and elevational directions. An observation in a given
The following section provides an overview of WSR- radar bin is an integrated returned power from all scat-
88D scan strategies and the associated beam geometry. ters within the sampling volume of the bin, with the
The objective analysis schemes mentioned above are scatters near the beam center weighted more than the
illustrated and contrasted in section 3. Strategies for scatters near the beam edges. At each range bin, the
combining data from multiple radars are presented in variance of the returned power is reduced to ⬃1 dB by
section 4, and a summary is provided in section 5. averaging a number of independent pulses (Sirmans
and Doviak 1973). The spatial density of radar data can
be discussed in terms of the distance between centers of
2. WSR-88D data resolution adjacent radar bins. In the radial direction the data
spacing of the WSR-88D reflectivity data is fixed at 1
The WSR-88Ds pulse electromagnetic energy into
km. In the azimuthal and elevational directions, how-
the atmosphere along a conical beam as the antenna
ever, the data spacing is a function of range (Fig. 2).
scans azimuthally in specified elevation angle steps.
The azimuthal data spacing is calculated by assuming
The resultant volume scan is in spherical (r, ␾, ␪e) co-
that radials are 1° apart and evenly distributed. There-
ordinates, where r is the slant range, ␾ is azimuth angle
fore, the distance between centers of adjacent radials is
from north, and ␪e is elevation angle from the horizon.
simply r公1 ⫺ cos(⌬␾), where ⌬␾ ⫽ 1°. ⌻he azimuthal
The electromagnetic power distribution in a beam fol-
data spacing is 1 km at a range of approximately 75 km
lows a Bessel function of second order (Doviak and
and increases linearly with range (Fig. 3a). The eleva-
Zrnic 1993):
tional data spacing is calculated in a similar way, but the

f 2共␣兲 ⫽ 再 8 J2关共␲D sin␣兲Ⲑ␭兴

关共␲D sin␣兲Ⲑ␭兴2
冎 2
, 共1兲
angle differences are those between two adjacent tilts
(prescribed by VCPs). Since adjacent tilt angles are not
uniform (Table 1), the elevational data spacing varies
where ␣ is the angular distance from the beam axis, J2 not only with range but also with elevation. The large
is the Bessel function of second order, D is the diameter gaps between higher tilts in VCP21 result in very poor
of the antenna reflector, and ␭ is the wavelength. vertical coverage (Fig. 3b). The data spacing is ⬎1 km
The operational WSR-88Ds scan in four different beyond a range of 30 km (15 km) above the fifth (sixth)

TABLE 1. Elevation angles (°) used in the four NWS operational VCPs for WSR-88Ds.
Tilt No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14
VCP 11 0.5 1.45 2.4 3.35 4.3 5.2 6.2 7.5 8.7 10.0 12.0 14.0 16.7 19.5
VCP 21 0.5 1.45 2.4 3.35 4.3 6.0 9.9 14.6 19.5
VCP 31 0.5 1.5 2.5 3.5 4.5
VCP 32 0.5 1.5 2.5 3.5 4.5
JANUARY 2005 ZHANG ET AL. 33

FIG. 2. WSR-88D radar power density function for (a) VCP 11

and (b) VCP 21. The brighter white areas indicate higher power
densities.

tilt. Significant data voids exist between higher tilts as

well as below the lowest tilt and above the highest tilt
(see black regions in Figs. 2a,b). These factors signifi-
cantly impact the spatial fidelity of the radar data, FIG. 3. WSR-88D data resolutions in azimuthal and elevational
directions for VCP 21.
which is apparent before and after the radar data have
been objectively analyzed to a grid.
300–400 km (Fig. 5a), yet the height of precipitating
3. Remapping radar data to Cartesian space echo is only approximately 7.5 km (Figs. 5b,c). A radar
bright band [resulting from melting hydrometeors
a. Case description (Atlas and Banks 1950)] is apparent at ⬃1.7 km above
radar, where reflectivity is ⬃45 dBZ at the center of the
Approximately 50 cases from different geographical
bright band and decreasing to ⬃30 dBZ within 250 m
regions and within each season were examined using
below and above (Figs. 5b,c).
the various analysis schemes. This section presents the
results from two representative cases, a convective
b. Radar bin volume mapping
storm event that occurred on 25 June 2002 in Indiana
(Fig. 4) and a wintertime stratiform precipitation event Four interpolation schemes were evaluated for
that occurred on 15 February 1998 in Arizona (Fig. 5). remapping radar data from their native spherical coor-
These two cases represent different weather regimes dinates to Cartesian coordinates. The first scheme, ra-
occurring on dissimilar spatial scales, as is evident in dar bin volume mapping (RBVM), simply fills in grid
their radar depictions. The convective storm case shows cells collocated within a given radar bin with the value
strong upright convective cells with horizontal scales on observed within the bin. A grid cell is determined to be
the order of 10–100 km with vertical scales up to 15 km within a radar bin if the center of the grid cell is con-
(Figs. 4b,c). The winter case, in contrast, shows much tained inside the volume of the radar bin. The volume
more horizontally homogeneous echoes (often termed of radar bins is assumed to be 1 km ⫻ 1° ⫻ 1°. 〈ll grid
“stratiform” rain) with large horizontal extents of cells that are not encompassed by any radar bin volume
34 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 22

FIG. 4. (a) Composite reflectivity and range–height indicator

(RHI) reflectivity plots along (b) 263° and (c) 122° azimuths. The
data are from KIWX (North Webster, IN) radar valid at 2036 FIG. 5. (a) Composite reflectivity and RHI reflectivity plots
UTC 25 Jun 2002. The thin white lines in (a) indicate constant along (b) 0° and (c) 270° azimuths. The data are from KIWA
ranges (every 50 km) and azimuths (every 45°). The bold white (Phoenix, AZ) radar valid at 0859 UTC 15 Feb 1998. The bold
lines in (a) indicate where the RHIs in (b) and (c) were obtained. white lines in (a) indicate where the RHIs in (b) and (c) were
obtained.
JANUARY 2005 ZHANG ET AL. 35

are flagged as missing (e.g., grid cells in the black re- Horizontal plots of reflectivity for the convective
gions in Fig. 2). Figures 4b, 4c, 5b, and 5c show vertical case provide a smooth depiction of the reflectivity field
cross sections of RBVM analyses from the two cases (Fig. 7). Nevertheless, arc-shaped discontinuities simi-
using a very fine Cartesian grid (50 m ⫻ 50 m ⫻ 10 m). lar to the winter brightband case can be present in the
Using this fine grid, the RBVM scheme depicts radar trailing stratiform region (Fig. 7). The black rings near
beam propagation and data distributions (assuming 4/3- the radar in both Figs. 6 and 7 indicate that gaps be-
effective earth radius model) and the radar sampling tween the higher tilts are related to the individual scan
limitations of the VCP 11 and VCP 21 scanning strat- strategy.
egies. A more complete depiction of convective reflec-
tivity structure is accomplished using VCP11 (Figs.
c. Nearest-neighbor mapping
4b,c) as compared to VCP 21 (Figs. 5b,c).
A horizontal cross section taken at 1.9 km above the The second scheme being examined, nearest neigh-
radar for the winter storm case shows ring-shaped ar- bor mapping (NNM), assigns the value of the closest
tifacts as a result of the relatively poor vertical sampling radar bin to grid cell, where distance is evaluated using
(Fig. 6). Radar bin sizes increase with range rapidly the location of the centers of the radar bins. The near-
from an approximate diameter of 0.3 km at a range of est-neighbor approach results in a horizontal cross sec-
25 km to 2.5 km at a range of 150 km (Fig. 5b). The tion (not shown) that is very similar to that produced by
1.9-km height intersects different rays at different dis- the RBVM method but with the characteristic of filling
tances. Sometimes grid points at the 1.9-km level are in data voids near the radar. Therefore, it suffers from
within radar bins that are strongly affected by the bright many of the same artifacts as the RBVM approach
band, while other times grid points are within radar (e.g., the ring-shaped artifacts). The ring artifacts of the
bins that are not strongly affected by the bright band. RBVM and NNM analysis schemes are not surprising
The high-reflectivity rings are associated with radar because these schemes are simply moving observation
bins that are centered near the bright band, while the values instead of estimating from trends in the data.
low-reflectivity rings are associated with radar bins that The problem is more pronounced in the vertical than in
are centered above or below the bright band. horizontal because the scales of atmospheric phenom-

FIG. 6. Horizontal cross section of reflectivity analysis at 1.9 km above the radar using the
RBVM method. The data are from KIWA and are valid at 0859 UTC 15 Feb 1998.
36 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 22

FIG. 7. Horizontal cross section of the reflectivity analysis at 4.7 km above the radar using
the RBVM method. The data are from KIWX and are valid at 2036 UTC 25 Jun 2002. The red
arrows indicate an arc-shaped discontinuity in the trailing stratiform region as a result of
RBVM.

ena are much smaller in the vertical than in the hori- Here w1 and w2 are the interpolation weights given to
zontal. the reflectivity observations below and above the grid
cell, respectively. The weights are determined by
d. Vertical interpolation
w2 ⫽ 共␪i ⫺ ␪o1 兲Ⲑ共␪o2 ⫺ ␪o1 兲 and 共6兲
The third scheme examined is a linear interpolation
scheme in the elevational direction combined with the
nearest neighbor scheme in the azimuthal and range w1 ⫽ 共␪o2 ⫺ ␪i兲Ⲑ共␪o2 ⫺ ␪o1 兲, 共7兲
directions. The vertical interpolation (VI) approach has
been used quite successfully in airborne Doppler radar where ␪i, ␪o1 , and ␪o2 represent elevation angles of the
applications (Jorgensen et al. 1996) where the antenna grid cell and the radar bins below and above, respec-
makes vertical scans. For small elevation angles (⬍20°), tively.
the elevational direction is approximately vertical By performing a linear interpolation in the eleva-
(hence the schemes name). The procedure for comput- tional direction, the vertical gradients are better pre-
ing the analysis value f ai at a given grid cell i is as served than with nearest-neighbor mapping. The verti-
follows. cal structure of convective storm cells is more coherent
than in the RBVM (Fig. 8a versus Fig. 4b) and with
1) Find the range, azimuth, and elevation at the center NNM analyses. Moreover, the ring-shaped discontinu-
of the grid cell i. ities in the brightband layer have been alleviated, ex-
2) Find two observations, f o1 and f o2 , on the two adja- cept at ranges corresponding to gaps between tilts in
cent tilts below and above the grid cell, respectively, VCP 21 (Figs. 9a,b versus Figs. 5b and 6), where the VI
and at the same range and azimuth as the grid cell. scheme relies on reflectivity observations above and
3) Compute the analysis value f ai using below the bright band to fill in the data voids. As a
result of the VI, low reflectivity rings appear in these
f ia ⫽ 共w1 f o1 ⫹ w2 f o2 兲Ⲑ共w1 ⫹ w2兲. 共5兲 gaps on the horizontal cross sections (Fig. 9b). To bet-
JANUARY 2005 ZHANG ET AL. 37

FIG. 8. (a) An RHI plot along 263° azimuth and (b) a horizontal
cross section at 4.7 km above the radar of the reflectivity analysis
FIG. 9. (a) An RHI plot along 0° azimuth and (b) a horizontal
using the VI scheme. The data are for the convective case. The red
cross section at 1.9 km above the radar of the reflectivity analysis
arrows in (b) indicate places where arc-shaped discontinuities re-
using the VI scheme. The data are for the winter stratiform case.
sulted from the RBVM approach but not from the VI.
The red arrows in (a) indicate places where reflectivity values
above and below the bright band are used to interpolate the re-
flectivity in the gaps.
ter preserve brightband layer structures, a fourth
scheme is employed and subsequently evaluated.
and w4, the interpolation weights given to the two ob-
e. Vertical and horizontal interpolation servations f o3 and f o4 , respectively, are determined by
The fourth scheme uses the VI scheme plus a hori- w4 ⫽ 共so3 ⫺ si兲/共so3 ⫺ so4 兲 and 共9兲
zontal interpolation between adjacent tilts that are
more than 1° apart. The analysis formula for the verti- w3 ⫽ 共si ⫺ so4 兲/共so3 ⫺ so4 兲. 共10兲
cal and horizontal interpolation (VHI) scheme is the
following: Here si represents horizontal distance between the ra-
dar (at the origin in Fig. 10) and the grid cell i; and so3
f ia ⫽ 共w1 f o1 ⫹ w2 f o2 ⫹ w3 f o3 ⫹ w4 f o4 兲/ and so4 represent horizontal distances between the radar
and data bins “3” and “4,” respectively (Fig. 10).
共w1 ⫹ w2 ⫹ w3 ⫹ w4兲. (8)
The horizontal interpolation, when applied in the
Here f o3 and f o4 represent the two reflectivity obser- data voids, recovered the brightband layer and further
vations along horizontal directions on the two adjacent alleviated the ring-shaped artifacts on horizontal cross
tilts below and above the grid cell (see Fig. 10); and w3 sections (Fig. 11). The 3D reflectivity analyses gener-
38 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 22

FIG. 10. An illustration of the VHI scheme. The grid point is

indicated by ⫹, and numbers 1–4 indicate the four radar bins
where the observations are used to compute the analysis value at
the grid point.

ated using the VHI approach also produce more real-

istic echo-top fields (Fig. 12). Ring artifacts often ap-
pear in echo-top fields derived from reflectivity in
spherical coordinates (Howard et al. 1997; Maddox et
al. 1999; Brown et al. 2000). They also exist in the re-
flectivity analysis using the RBVM (Fig. 12a) and the
NNM (not shown) approaches. The VHI approach,
which makes uses of trends (with a linear approxima-
tion) in the data, is able to alleviate the ring-shaped
discontinuities significantly (Fig. 12b).
It was discovered that the VHI scheme is not as well
suited for the convective echoes characterized by sig-
nificant variability in the horizontal. As shown in Fig.
13, for example, the VHI approach produces unwanted FIG. 11. Same as in Fig. 9 except for the VHI scheme. The red
artifacts near the edges of the upright convective cores. arrows in (a) indicate places where artifacts due to the vertical
part of the interpolation scheme remain but to a much lesser
Based on the experiments above and the detailed ex- degree than in the VI scheme.
amination of several cases (convective and stratiform),
the VI approach is suitable for most situations. The
exception is the occurrence of a bright band, where the mapped onto a common Cartesian grid. The remapped
VHI approach preserves more of the continuous nature reflectivity fields from multiple radars are then com-
of the brightband layer and avoids ring-shaped artifacts bined, or “mosaicked,” to produce a unified 3D reflec-
in the analysis. Therefore, the objective analysis ap- tivity grid. There are many areas across the United
proach requires adaptation to the spatial scales (hori- States, especially at mid- to upper levels of the tropo-
zontal and vertical) of the observed weather. Objective sphere, where multiple radars overlap in coverage
brightband identification schemes such as the one de- (Maddox et al. 2002). The reflectivity value for each
veloped by Gourley and Calvert (2003) can be used to grid cell in a mosaicked grid i is obtained by
determine the existence of a bright band and thus de-
termine a proper objective analysis approach in real time. f m共i兲 ⫽ 兺
n⫽1,Nrad
␣n f na共i兲 冒 兺
n⫽1,Nrad
␣n. 共11兲

4. Mosaic schemes Here f m(i) is the mosaicked reflectivity value at the

grid cell, i; f an(i) is an analysis value at the grid cell from
a. Multiradar mosaic the nth radar; ␣n is a weight given to the analysis value
Through utilization of one or more of the objective f an(i); and Nrad represents the total number of radars
analysis approaches discussed in the previous section, that have an analysis value at the grid cell. If Nrad ⫽ 0,
the reflectivity fields from individual radars are re- then the grid cell is not covered by any radar and a
JANUARY 2005 ZHANG ET AL. 39

FIG. 13. Same as in Fig. 9a except for the VHI scheme. Note that
artifacts (indicated by the red circle) are introduced by the hori-
zontal interpolation.

in the mosaicked field. Figures 14a and 14b show com-

posite reflectivity fields for two adjacent radars (KIWX
at North Webster, Indiana, and KLOT at Romeoville,
Illinois) that share a common area and were valid near
the same time. The KIWX composite reflectivity values
are significantly higher than those of KLOT in the ma-
jority of the domain shared by both radars (Figs. 14a,b).
The mean difference (KIWX ⫺ KLOT) computed at
the equidistant grid cells along the red line in Fig. 14
between the two fields is ⫹7 dBZ. Gourley et al. (2003)
has shown that calibration differences among WSR-
88Ds often exceed 5 dBZ. The calibration differences
result in discontinuities in the nearest-neighbor mosaic
along the equidistant lines, making the nearest-
neighbor scheme suboptimal for mosaicking reflectivity
from multiple radars. Discontinuities could also arise
from other causes, such as the presence of a bright band
that is sampled differently by the two radars.
The second mosaicking method simply uses the maxi-
FIG. 12. Echo-top fields (in km above radar level) derived from
the 3D reflectivity analyses (KIWA, 0859 UTC 15 Feb 1998) using
mum reflectivity value among the multiple observations
the (a) RBVM and (b) VHI approaches. that cover the same grid cell. This method does not
involve smoothing and it retains the highest reflectivity
intensities in the data fields. The method possibly miti-
missing flag is assigned to f m(i). If Nrad ⫽ 1, then the gates attenuation losses due to intervening regions of
grid cell is covered by one radar only and the analysis intense rainfall. Therefore, it is better suited for creat-
value from that radar is assigned to the grid cell [i.e., ing derivative products such as composite reflectivity
f m(i)⫽f a1(i)]. If Nrad ⬎ 1, then a weighted mean of fields. However, this method would provide a biased
multiple-radar analysis values is assigned to the grid estimate toward radars that provide higher reflectivity
cell. Several weighting methods are considered and values as a result of calibration differences.
tested to determine an optimal weighting methodology The third mosaic method considered is a weighted
for combing radar data. mean, whereby the weight is based on the distance be-
The first mosaic weighting method is the “nearest tween an individual grid cell and the radar location.
neighbor,” in which the analysis value from the closest Two weighting functions were tested; both monotoni-
radar is assigned to the grid cell. The nearest-neighbor cally decrease with range (Fig. 15). The weighting func-
method does not impose any smoothing when creating tion reflects the confidence level that the radar obser-
a mosaic from multiple radars. However, discontinu- vation is representative of the convective-scale storm
ities may appear at the equidistant lines between radars structure. Due to the beam spreading, the size of the
40 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 22

FIG. 14. Composite reflectivity fields observed by the (a) KLOT and (b) KIWX radars for the same region at 2036
UTC 25 Jun 2002 and the mosaicked composite reflectivity fields using the (c) flat and (d) steep mosaic weighting
functions shown in Fig. 15. The red lines in (a) and (b) indicate the equidistant line between the two radars.

radar resolution volume increases squarely with range. observations from distant radars. Mosaicked composite
Subsequently, the reflectivity fields have less finescale reflectivity fields from KIWX and KLOT radars were
structure at far ranges than at ranges closer to the radar. derived using the two weighting functions (Figs. 14c,d).
Therefore, the observations at far range should receive The mosaicked composite reflectivity field using the
smaller weights than those from near the radar when flat weighting function scheme (Fig. 14c) showed much
mosaicking. The first weighting function has a charac- lower reflectivities than those depicted in the KIWX
teristically steep shape (dashed line in Fig. 15) where observations (Fig. 14b), especially in regions near the
the weight decreases with range very rapidly. In con- KIWX radar. The lesser reflectivity intensities in this
trast, the second weighting function (solid line in Fig. analysis are due to the fact that lower reflectivity values
15) is relatively flat, so that more weight is given to from the KLOT radar received significant weights dur-
JANUARY 2005 ZHANG ET AL. 41

5. Summary

Four interpolation approaches and three mosaic

methods are evaluated in this paper. These techniques
are candidates for converting radar reflectivity data
from multiple sources in their native (spherical) coor-
dinates onto a multiradar, mosaicked 3D Cartesian
grid. Approximately 50 cases from different geographi-
cal regions and different seasons were evaluated and
examples of two common weather regimes (spring/
summer convective storms and wintertime stratiform
precipitation systems) are presented. Based on exten-
sive case studies, it was found that the vertical interpo-
lation scheme provides the most physically realistic mo-
saic for convective storms. For stratiform precipitation
with the presence of a bright band, an additional hori-
zontal interpolation scheme is required to reconstruct
FIG. 15. Two examples of mosaic weighting functions. the horizontally extended layers in the data voids, es-
pecially between upper tilts, as in VCP 21. When com-
bining multiple-radar reflectivity fields onto a single 3D
ing the mosaicking analysis. In addition, due to the
grid, a distance-weighted mean scheme is preferred to
beam spreading, reflectivity observations at far ranges
produce a continuous and representative field. The
contain features with coarser scale than reflectivity ob-
shape of the weighting function is important for retain-
servations close to the radar. If a flat weighting function
ing the spatial reflectivity gradients and intensities.
is applied, the far-range observations can dampen the
By using the vertical interpolation and vertical/
finescale gradients associated with near-range observa-
horizontal interpolation approaches together with a dis-
tions. Conversely, the steep weighting function retains
tance-weighted mosaicking scheme, multiple-radar re-
more details on reflectivity structures by allowing more
flectivity fields can be integrated onto a common 3D
contributions from near-range observations. The latter
Cartesian grid. These mosaics contain high-resolution
scheme is thus preferred in this case for mosaicking
features of the raw data, but in a manageable size with
reflectivity data from multiple radars due to its capa-
the advantage of a Cartesian coordinate system. The
bility of representing finescale storm structures while
high-resolution mosaic grid provides useful datasets for
preserving their magnitudes.
further applications, such as quantitative precipitation
estimation, severe storm detection algorithms, and data
b. Data voids
assimilation for convective-scale numerical weather
Figure 2 shows the existence of gaps between the predictions.
higher tilts as well as the data voids above the highest The VI and VHI schemes have the advantage of low
beam (i.e., the “cone of silence”) and below the lowest computational costs. However, the resultant reflectivity
beam. All the spherical-to-Cartesian remapping schemes analyses are not very suited for applications that re-
discussed in section 4 do not extrapolate reflectivity quire high-order derivatives of the analyses. This is due
values into the data void regions below the beam cov- to the fact that linear interpolations used in the scheme
erage shown in Fig. 2. While multiple-radar mosaic can can introduce discontinuities in derivative fields of the
fill in the upper part of one radar’s “cone of silence” analyses. A smoothing or filtering of the mosaic grid
using observations from other radars nearby, significant data is recommended before such applications. With
data voids still exist in the lower altitudes (Maddox et the Cartesian grid, filtering of the data is straightfor-
al. 2002). Without filling in the data voids, the applica- ward in comparison to spherical coordinates.
tion of the 3D mosaic grid to numerical weather models Even with the interpolations and multiple-radar mo-
is limited since the information near the surface is criti- saic, data voids still exist such as with the cone of silence
cal to model initialization. An additional radar net- above the radar and in regions below the lowest beam.
work, such as the FAA Terminal Doppler Weather Ra- To reduce the data voids, the FAA TDWR data will be
dar (TDWR) network, can potentially reduce the data integrated into the mosaic for better lower-atmosphere
voids in areas near airports and large metropolitan ar- coverage, especially over metropolitan areas. New stud-
eas (Vasiloff 2001). Vertical profiles of reflectivity ies are under way to assess the potential of using ver-
(VPRs; e.g., Joss and Waldvogel 1990; Berne et al. tical profiles of reflectivity to extrapolate reflectivity
2004) can be derived from the observations at the close observations at far ranges down to the earth’s surface.
ranges. The profiles can be used to approximate verti- Future research also includes development of a syn-
cal reflectivity structures at far ranges down to the chronization scheme for handling time and sampling
earth’s surface. differences among multiple radars’ observations. These
42 JOURNAL OF ATMOSPHERIC AND OCEANIC TECHNOLOGY VOLUME 22

studies will enhance the 3D mosaic grid and increase its monitoring thunderstorm life cycles. Wea. Forecasting, 12,
usefulness for forecasters and for the inclusion of high- 166–174.
Jorgensen, D. P., P. H. Hildebrand, and C. L. Frush, 1983: Feasi-
resolution 3D mosaic fields in numerical models. bility test of an airborne pulse-Doppler meteorological radar.
J. Appl. Meteor., 22, 744–757.
Acknowledgments. The authors would like to thank ——, T. Matejka, and J. D. DuGranrut, 1996: Multi-beam tech-
Dr. Dave Jorgensen and Dr. Robert Maddox and the niques for deriving wind fields from airborne Doppler radars.
J. Meteor. Atmos. Phys., 59, 83–104.
reviewers for their comments that helped immensely in
Joss, J., and A. Waldvogel, 1990: Precipitation measurement and
the preparation of this manuscript. Major funding for hydrology. Radar in Meteorology, D. Atlas, Ed., Amer. Me-
this research was provided under the Aviation Weather teor. Soc., 577–606.
Research Program NEXRAD Algorithms Product De- Maddox, R. A., D. S. Zaras, P. L. Mackeen, J. J. Gourley, R.
velopment Team (NAPDT) MOU, and partial funding Rabin, and K. W. Howard, 1999: Echo height measurements
with the WSR-88D: Use of data from one versus two radars.
was provided under NOAA-OU Cooperative Agree- Wea. Forecasting, 14, 455–460.
ment NA17RJ1227. ——, J. Zhang, J. J. Gourley, and K. W. Howard, 2002: Weather
radar coverage over the contiguous United States. Wea. Fore-
casting, 17, 927–934.
REFERENCES Miller, L. J., C. G. Mohr, and A. J. Weinheimer, 1986: The simple
rectification to Cartesian space of folded radial velocities
Askelson, M. A., J.-P. Aubagnac, and J. M. Straka, 2000: An from Doppler radar sampling. J. Atmos. Oceanic Technol., 3,
adaptation of the Barnes filter applied to the objective analy- 162–174.
sis of radar data. Mon. Wea. Rev., 128, 3050–3082. Mohr, C. G., and R. L. Vaughn, 1979: An economical procedure
Atlas, D., and H. C. Banks, 1950: A virtual echo-layer above the for Cartesian interpolation and display of reflectivity factor
bright band. J. Atmos. Sci., 7, 402–403. data in three-dimensional space. J. Appl. Meteor., 18, 661–
Berne, A., G. Delrieu, H. Andrieu, and J.-D. Creutin, 2004: In- 670.
fluence of the vertical profile of reflectivity on radar- OFCM, 2003: Doppler radar meteorological observations. Part A:
estimated rain rates at short time steps. J. Hydrometeor., 5, System concepts, responsibilities, and procedures. Federal
296–310. Meteorological Handbook 11, Office of the Federal Coordi-
Bluestein, H. B., and S. G. Gaddy, 2001: Airborne pseudo-dual- nator for Meteorological Services and Supporting Research,
Doppler analysis of a rear-inflow jet and deep convergence 34 pp.
zone within a supercell. Mon. Wea. Rev., 129, 2270–2289. Rinehart, R. E., 2001: Radar for Meteorologists. 3d ed. Rinehart,
Brown, R. A., V. T. Wood, and D. Sirmans, 2000: Improved WSR- 428 pp.
88D scanning strategies for convective storms. Wea. Forecast- Shapiro, A., P. Robinson, J. Wurman, and J. Gao, 2003: Single-
ing, 15, 208–220. Doppler velocity retrieval with rapid-scan radar data. J. At-
Doviak, R. J., and D. S. Zrnic, 1993: Doppler Radar and Weather mos. Oceanic Technol., 20, 1758–1775.
Observations. 2d ed. Academic Press, 562 pp. Sirmans, D., and R. J. Doviak, 1973: Doppler radar with polarized
Droegemeier, K. K., and Coauthors, 2002: Project CRAFT: A test diversity. J. Atmos. Sci., 30, 162–174.
bed for demonstrating the real time acquisition and archival Trapp, R. J., and C. A. Doswell, 2000: Radar data objective analy-
of WSR-88D Level II Data. Preprints, 18th Int. Conf. on sis. J. Atmos. Oceanic Technol., 17, 105–120.
Interactive Information Processing Systems (IIPS) for Meteo- Vasiloff, S., 2001: Improving tornado warnings with the Federal
rology, Oceanography, and Hydrology, Orlando, FL, Amer. Aviation Administration’s Terminal Doppler Weather Ra-
Meteor. Soc., 136–139. dar. Bull. Amer. Meteor. Soc., 82, 861–874.
Fulton, R., 1998: WSR-88D polar to HRAP mapping. Tech. Weygandt, S. S., A. Shapiro, and K. K. Droegemeier, 2002a: Re-
Memo., Hydrology Research Laboratory, Office of Hydrol- trieval of model initial fields from single-Doppler observa-
tions of a supercell thunderstorm. Part I: Single-Doppler ve-
ogy, National Weather Service, Silver Spring, MD, 33 pp.
locity retrieval. Mon. Wea. Rev., 130, 433–453.
Gao, J., M. Xue, A. Shapiro, and K. K. Droegemeier, 1999: A
——, ——, and ——, 2002b: Retrieval of model initial fields from
variational method for the analysis of three-dimensional wind
single-Doppler observations of a supercell thunderstorm.
fields from two Doppler radars. Mon. Wea. Rev., 127, 2128–
Part II: Thermodynamic retrieval and numerical prediction.
2142. Mon. Wea. Rev., 130, 454–476.
Gourley, J. J., and C. M. Calvert, 2003: Automated detection of Wurman, J., and S. Gill, 2000: Finescale radar observations of the
the bright band using WSR-88D radar data. Wea. Forecast- Dimmitt, Texas (2 June 1995), tornado. Mon. Wea. Rev., 128,
ing, 18, 585–599. 2135–2164.
——, B. Kaney, and R. A. Maddox, 2003: Evaluating the calibra- Zhang, J., 1999: Moisture and diabatic initializations based on
tions of radars: A software approach. Preprints, 31st Int. radar and satellite observations. Ph.D. dissertation, Univer-
Conf. on Radar Meteorology, Seattle, WA, Amer. Meteor. sity of Oklahoma, 194 pp.
Soc., 459–462. ——, and T. Gal-Chen, 1996: Single-Doppler wind retrieval in the
Heymsfield, G. M., 1976: Statistical objective analysis of dual- moving frame of reference. J. Atmos. Sci., 53, 2609–2623.
Doppler radar data from a tornadic storm. J. Appl. Meteor., Ziegler, C. L., E. N. Rasmussen, T. R. Shepherd, A. I. Watson,
15, 59–68. and J. M. Straka, 2001: The evolution of low-level rotation in
Howard, K. W., J. J. Gourley, and R. A. Maddox, 1997: Uncer- the 29 May 1994 Newcastle–Graham, Texas, storm complex
tainties in WSR-88D measurements and their impacts on during VORTEX. Mon. Wea. Rev., 129, 1339–1368.