Three-Dimensional Cloud Volume Reconstruction from the Multi-angle Imaging SpectroRadiometer

Abstract

Characterizing 3-D structure of clouds is needed for a more complete understanding of the Earth’s radiative and latent heat fluxes. Here we develop and explore a ray casting algorithm applied to data from the Multi-angle Imaging SpectroRadiometer (MISR) onboard the Terra satellite, in order to reconstruct 3-D cloud volumes of observed clouds. The ray casting algorithm is first applied to geometrically simple synthetic clouds to show that, under the assumption of perfect, clear-conservative cloud masks, the reconstruction method yields overestimation in the volume whose magnitude depends on the cloud geometry and the resolution of the reconstruction grid relative to the image pixel resolution. The method is then applied to two hand-picked MISR scenes, fully accounting for MISR’s viewing geometry for reconstructions over the Earth’s ellipsoidal surface. The MISR Radiometric Camera-by-camera Cloud Mask (RCCM) at 1.1-km resolution and the custom cloud mask at 275-m resolution independently derived from MISR’s red, green, and blue channels are used as input cloud masks. A wind correction method, termed cloud spreading, is applied to the cloud masks to offset potential cloud movements over short time intervals between the camera views of a scene. The MISR cloud-top height product is used as a constraint to reduce the overestimation at the cloud top. The results for the two selected scenes show that the wind correction using the cloud spreading method increases the reconstructed volume up to 4.7 times greater than without the wind correction, and that the reconstructed volume generated from the RCCM is up to 3.5 times greater than that from the higher-resolution custom cloud mask. Recommendations for improving the presented cloud volume reconstructions, as well as possible future passive remote sensing satellite missions, are discussed.

1. Introduction

Characterizing 3-D cloud structures, such as shape and volume, is critical in understanding the distribution of radiative and latent heat fluxes on Earth (e.g., References [1,2]). Compared to traditional 1-D or 2-D treatments, 3-D cloud structures have a distinguishing impact on not only the forward calculations of the fluxes (e.g., References [3,4,5,6]), but also the inverse problem of retrieving cloud optical and microphysical properties from radiation measurements (e.g., References [7,8,9,10]). Indeed, their impact is evident in the spectral, textural, and angular distributions of the outgoing radiation field observed by satellites (e.g., References [11,12]). Among many cloud properties to be characterized in 3-D, cloud shape and volume could be especially useful in setting 3-D cloud boundaries for radiance calculations, such as for cloud radiative effects and 3-D scattering.

While our Earth-observing satellite system has advanced over the past several decades, our ability to retrieve 3-D cloud structures is incomplete, and there exist only a limited number of studies that go beyond 2-D cloud structural properties. Zinner et al. (2006) [13] showed the reconstruction of 3-D stratocumulus structures from high-spatial-resolution radiance fields observed by the Compact Airborne Spectrographic Imager (CASI, 15-m resolution). While promising, its application was limited to small solar zenith angles (no shadows allowed) and single-layered stratiform clouds, and the method has not yet found applications using satellite data. In a study by Barker et al. (2011) [14], an algorithm was developed and assessed for reconstructing the 3-D distribution of clouds from passive satellite imagery and collocated 2-D nadir profiles of cloud properties inferred synergistically from lidar, cloud radar, and imager data. It has found applications to solar radiative transfer calculations with the A-Train [15,16], but has limitations for optically thick cloud fields and relies on a statistical approach for reconstructions rather than direct measurements. Ground-based scanning cloud radars are capable of retrieving 3-D cloud structures [17,18], but limitations, due to their immobility relative to satellites and radar beam attenuation, will be at play in covering wider areas, such as oceans, and reconstructing vertically developed clouds, respectively.

Bearing in mind the current status of the research area, we present a method of reconstructing 3-D cloud volumes using data from the Multi-angle Imaging SpectroRadiometer (MISR) [19]. MISR, launched in 1999 onboard NASA’s Terra satellite, views the Earth from a sun-synchronous orbit at nine widely spaced view angles. The camera designations (viewing zenith angles) in the along-track direction are from the most forward view DF (${70.5}^{°}$) to CF (${60.0}^{°}$), BF (${45.6}^{°}$), AF (26.1${}^{°}$), the nadir view AN ($0^{°}$), AA ($-{26.1}^{°}$), BA ($-{45.6}^{°}$), CA ($-{60.0}^{°}$), and then to the most aft-ward (backward) view DA ($-{70.5}^{°}$). MISR takes about seven minutes to look at the same scene on the Earth’s surface from the DF camera to the DA camera. MISR has four spectral channels representing blue (446 nm), green (558 nm), red (672 nm), and near infrared (866 nm). The measured radiances are stored with the highest spatial resolution of 275 m per pixel on the Space-oblique Mercator (SOM) projection grid, with a cross-track swath of around 400 km. For geo-registration and storage, MISR data are broken into a series of predefined, uniformly-sized boxes along the ground track called blocks. Each orbital path is divided into 180 blocks measuring 563.2 km (cross-track) x 140.8 km (along-track). MISR has global coverage with 233 orbits in the period of sixteen days. Its current operational period of 18 years and ongoing, together with its multi-angle viewing capability, presents a unique opportunity to study the 3-D distribution of clouds around the globe from a daytime orbit over a substantial period of time.

There have been a few studies that explore the use of MISR for retrieving cloud geometry. Kassianov et al. (2003) [20] generated their own cloud masks and retrieved the geometrical thickness of cloud from AN image pixels’ radiances using simple models, which was then “forced” to agree with MISR’s other multi-angle observations. Seiz and Davies (2006) [21] and Cornet and Davies (2008) [22] presented 2-D cross-sectional reconstructions using stereo-photogrammetric methods. However, reconstructions of 3-D cloud shape and volume from MISR observations have yet to be explored. Here, we present a ray casting method for reconstructing 3-D volumes of clouds from MISR data.

Section 2 describes the general principle behind the reconstruction method that uses a ray casting algorithm. Section 3 contains assessment of the method via simulations using synthetic clouds. Section 4 contains list of products and methods for reconstruction from MISR data. Section 5 presents cloud volume reconstruction results for two hand-picked MISR scenes. Section 6 discusses limitations and implications of the results and concludes the paper.

2. General Principle behind the Reconstruction Method

The reconstruction method aims at 3-D cloud volume reconstruction from 2-D satellite imagery. Focusing on the use of MISR, it employs a ray casting algorithm that rewinds the process of MISR’s atmospheric observations and data storage. Through a pre-determined 3-D reconstruction domain consisting of voxels (minimum computational units), rays are cast from the geo-registered locations of the cloudy pixels on the surface to the sensor’s location according to the sensor’s viewing zenith and azimuth angles, hence representing the radiance paths from the observed cloud to the satellite. Then, a series of calculations with the geometric properties of the voxels and the rays determine the space overlapped at all view angles, which is deemed to be the maximum space within the domain that the cloud occupies. Figure 1 gives a schematic diagram of this process.

Due to the limited number of views available in practice, this reconstruction method is expected to always overestimate the cloud volume (assuming that the cloud is perfectly detected in each image), especially at the cloud top and bottom, as evident by comparing the red hexagon with the cloud figure in Figure 1. To reduce excessive overestimation at cloud top, the standard cloud-top height product derived via a stereo-photogrammetric technique by the MISR science team [23,24,25] is considered as a constraint. Overestimation at cloud bottom could be reduced by including the lifting condensation level from reanalysis data or ground-based observations of cloud-base height. However, we do not include the application of the cloud-base height here, strictly focusing on the use of MISR data only.

The ray casting method presumes that any significant parcel of cloud in the atmosphere, unless occluded by another parcel, must be visible at all view angles. In other words, any parcel of cloud that does not appear in all view angles is either occluded or has negligible size and optical depth for the purpose of volume reconstruction. This assumption helps exclude pixels that are falsely determined cloudy in the cloud masks due to sun glint, which predominantly appears in only one or a few MISR view angles [26].

3. Simulations with Synthetic Clouds

To understand and quantify the errors and uncertainties expected with the reconstruction method when applied to MISR data, simulations were run with synthetically generated clouds. Synthetic satellite images were derived from these clouds through the process of image formation, and then were used as inputs into the cloud volume reconstruction simulations. The reconstructed volumes were then compared with the original synthetic-truth clouds.

Image formation was conducted as follows. First, from a 3-D synthetic cloud field, a set of 2-D binary cloud masks was derived in a way that imitates MISR’s retrievals with its nine along-track view angles of ±70.5°, ±60.0°, ±45.6°, ±26.1°, and 0° (nadir), the orbital height of 705 km, and the pixel resolution of 275 m. The ground instantaneous field of view (GIFOV) was set to be the square area occupied by the pixel, and the ground sampling interval (GSI) was set to be the same as the pixel resolution. The geo-registered locations of the pixels were set to be the same in all view angles, as in MISR data. The satellite camera location for each pixel at each view angle was set to imitate MISR’s push-broom technique, such that in the along-track axis the camera’s relative location with respect to each pixel’s geolocation remained the same according to the view angle, while in the cross-track axis its location remained fixed at the image’s center. This gave rise to changes in view angle across the swath. The dimensions of the images were set such that they included the clear ends of the cloud. For ray casting, the method of voxel boundary intersections was used: If the ray, cast from the geo-registered location of a pixel to the corresponding satellite camera location, passed through the boundary of any cloudy voxel, the pixel was determined cloudy. The voxel boundaries were determined in terms of geometric planes corresponding to the sides of the voxel cuboids, either perpendicular or parallel to the plane ground. The ray casting was repeated for all geo-registered cloudy pixels at all nine view angles to produce perfect synthetic cloud mask images.

Reconstructions were run with the cloud masks generated from the image formation as inputs. Under the same geometric parameters used in image formation, rays were cast from the geolocations of the cloudy pixels through the reconstruction domain to the corresponding satellite locations, and every voxel whose boundary was intersected by a ray was determined to be cloudy (cf. Figure 1). This was repeated for each of the nine view angles, and then only the voxels that had been determined cloudy at all angles were eventually determined cloudy in the final reconstructed volume.

To simplify the process, all simulations were run with a geometrically straight orbital path as opposed to MISR’s elliptical path. They were also conducted under the plane ground assumption, where the pixels were geo-registered to a ground-level plane instead of to terrain or an ellipsoidal surface as in MISR data. While there are many other different factors that affect the accuracy of the reconstructions, such as the number of view angles, the pixel resolution in which the satellite observation is made, and elements of cloud detection (instrument sensitivities, cloud optical depths, etc.), the synthetic clouds were suitable in assessing the reconstruction method in two major factors: The geometry of observed clouds and the resolution of reconstruction voxels relative to the image pixel resolution.

3.1. Cloud Geometry

To examine the behavior of the reconstructions with respect to cloud shape, simulations were run on four types of geometrically simple synthetic clouds: Cubic cloud, cross-track row cuboid cloud, along-track row cuboid cloud, and vertical column cuboid cloud. These basic cloud shapes were chosen to represent a unit parcel and the three fundamental axes (cross-track, along-track, and altitude) in three-dimensional geometry. Understanding the reconstructions from these basic shapes will provide a framework for understanding reconstructions from more complicated shapes, such as anvil clouds and distinctly layered clouds.

The synthetic clouds were put on the center of the image domain during the image formation, so that the nadir line of sight from the satellite aligns with the center of the cloud volume. The cloud-base height was set as 1650 m, within the usual altitude range (1–3 km) of real-world cumulus clouds. The cloud voxel resolution was set as 1650 m (i.e., the smallest unit volume of a reconstructed cloud was 1650³ m³). This number was chosen to ensure that the cloud is larger than the pixel resolution of 275 m and preclude excessive overestimation of cloudiness by otherwise too tiny (sub-pixel) clouds in the image formation of the perfect cloud masks [27]. The reconstruction voxel resolution was set as 275 m, same as the pixel resolution.

Figure 2 shows 3-D views of an example comparison between original and reconstructed cloud volume. Figure 3 shows the overall results for different shapes and lengths of the original cloud, in terms of the factor of overestimation (i.e., the ratio of reconstructed volume to original volume) as a function of the length of the longest side of the original cloud. As expected in our discussion of Figure 1, both Figures show that in all cases the reconstruction yields overestimation of cloud volume. Overestimation is particularly large in the cases of the along-track row cuboid cloud, resulting from the greater building of prisms growing with the cloud length at cloud top and bottom. The factor is around 2 for the cubic cloud, increases somewhat linearly at a rate of 0.1 km⁻¹ as the cloud length increases in the along-track direction, and decreases asymptotically to around 1.5 or 1.3 as the cloud length increases in the cross-track or the vertical direction, respectively. While the numerical values are particular to these simulation cases, it is possible to infer that, in general, the less horizontally spread out in the along-track direction and the more vertically developed the cloud, the less the overestimation in the cloud volume. Hence, this reconstruction method would yield the best results when applied to regions with vertically well-developed and horizontally narrow cumulus clouds, particularly if they were aligned in the cross-track direction. Note that having more oblique views beyond the MISR view angles would lead to less of the prism building at cloud top and bottom, thereby giving better reconstructions.

3.2. Reconstruction Voxel Resolution relative to the Image Pixel Resolution

To assess reconstructions in light of the reconstruction voxel resolution relative to the image pixel resolution, simulations were run for different values and the constant pixel resolution of 275 m. Cubic clouds with different sizes were tested. All other parameters, such as the cloud-base height and the viewing geometry, remained the same as in the previous simulations.

Figure 4 shows 3-D views of an example set of results for the cubic cloud with an edge = 1650 m. Figure 5 shows the overall results for cubic clouds with different sizes in terms of the factor of overestimation in cloud volume. Both Figures show that overestimation of cloud volume tends to decrease as voxel resolution becomes smaller. In other words, the smaller the reconstruction voxel resolution in comparison to the pixel resolution, the smaller the overestimation. The factor of overestimation never reaches 1, however, due to the method’s limitation discussed in Section 2. For most of the cubic cloud sizes tested, overestimation is more or less stable when the reconstruction voxel resolution is smaller than or equal to the pixel resolution. This is a useful reference for determining the reconstruction voxel resolution in the reconstructions from MISR data.

Figure 5 also shows that the greater the size of the original cloud with respect to the image pixel resolution, the smaller the overestimation of cloud volume. The relatively extensive overestimations at the smaller sizes of the original cloud (as shown especially by the cubic cloud with edge = 275 m) are due to the clear-conservative nature of the image formation, wherein smaller clouds still contribute the same to the cloudy determination of larger pixels (i.e., a pixel that is only partially cloudy is counted as fully cloudy) [27]. The overall bumpiness (lack of smoothness) in the curves is due to mismatching of resolutions, known as the resolution effect [27]. The 2-D resolution effect occurs between the original cloud and the synthetic satellite image pixels, and the 3-D resolution effect occurs between the synthetic images and the reconstruction voxels. In summary, we see that the reconstruction method works best when the observed cloud is large and the reconstruction voxel resolution is small in comparison to the pixel resolution.

4. Reconstruction from MISR Data: Methodology

Cloud volume reconstructions from MISR data entail details that are different from the synthetic simulations described in Section 3. In particular, (1) MISR data are geo-registered to a non-flat surface that accommodates the Earth’s surface curvature; (2) cloud mask classifications for cloudy and clear pixels may not be perfect; and (3) possible cloud movement due to wind and cloud morphing over the short time intervals between view angles call for the need for wind correction. These issues are addressed in the following subsections.

4.1. Reconstruction Domain

The reconstruction domain of voxels is defined in a way that reflects the geometry behind the geo-registration of MISR data. MISR uses two map projections for its Level 1 and Level 2 products: Ellipsoidal and terrain [28]. In this study, only the data in the ellipsoidal projection were used. The MISR ellipsoidal projection uses the World Geodetic System 1984 (WGS84) ellipsoid, with the equatorial radius of 6378137 m, the flattening factor of 298.257223563, and consequently the polar radius of 6356752 m [29]. The data are registered onto a rectangular grid, with the axes named SOM X and SOM Y, pre-defined by the space-oblique Mercator (SOM) projection of the surface of the ellipsoid. Accordingly, we defined the reconstruction domain in terms of the same SOM X and SOM Y as the horizontal axes, and the altitude as the vertical axis. Each voxel therefore has plane sides and ellipsoidal top and bottom. Hence, it forms a shape that is not strictly a rectangular cuboid but one distorted upward whose volume increases as the altitude increases (Figure 6). This is similar to how 3-D grid cells are defined in climate models. In this way, the curvature of the Earth’s ellipsoidal surface was reflected in the reconstruction.

Ellipsoidal reconstruction is employed as opposed to the plane-ground assumption used in Section 3, because the difference could lead to a noticeable difference in the reconstructed volumes. It was calculated, for the reconstruction voxel resolution = 1.1 km, that the ellipsoidal reconstruction domain, in comparison to the plane-ground reconstruction domain, has a voxel grid distortion reaching around 10 m and the block grid distortion reaching around 870 m at the altitude of 20 km. Here, the voxel grid distortion is defined as the difference between the reconstruction voxel resolution (1.1 km) and the straight-line distance from a voxel center to another adjacent at 20-km altitude. The block grid distortion is defined as the difference between the MISR block’s size in SOM Y (563.2 km) and the straight-line distance at the same altitude from one edge to another of the ellipsoidal reconstruction domain over the same block.

4.2. MISR Science Data

Multiple MISR science data products generated by the MISR science team and available on NASA’s EOSDIS website were used in the reconstructions. Their names and version numbers are listed as follows:

• MISR Radiometric Camera-by-camera Cloud Mask (RCCM) V004
• MISR Level 1B2 Ellipsoid Data V003
• MISR Ancillary Geographic Product (AGP) V001
• MISR Geometric Parameters (GMP) V002
• MISR Level 2 TOA/Cloud Height and Motion Parameters (TC_CLOUD) V001

The RCCM was one of the two input cloud masks used in the cloud volume reconstruction. The algorithm and theoretical basis for the derivation of the RCCM are described in the MISR Algorithm Theoretical Basis Document (ATBD) for Level 1 Cloud Detection [30], with the update for over ocean given by Zhao and Di Girolamo (2004) [26]. As its name suggests, the RCCM is independently derived for each of the nine MISR cameras. It has a pixel resolution of 1.1 km and 5 different levels of classification: High-confidence cloud, low-confidence cloud, low-confidence clear, high-confidence clear, and no retrieval. In this study, the RCCM was converted to a binary cloud mask depending on inputs for the classification levels. The RCCM pixels were assumed to have GIFOV and GSI equal to 1.1 km. Although in reality the GIFOV of the RCCM increases at oblique angles up to 1.53 km in the along-track direction, the effect of cloud overestimation by the oblique extension is negligible [26].

The L1B2 Ellipsoid Data provides the measured radiances for MISR’s red, green, and blue (RGB), and near infrared spectral channels. The RGB images were used in selecting scenes appropriate for reconstruction and the derivation of the corresponding custom cloud mask at 275-m resolution, which were another input cloud mask to the cloud volume reconstruction. The Ellipsoid Data also provides the reference time for each block at each view angle, which was used in the cloud mask wind correction (Section 4.4).

The AGP provides the geographical locations for each pixel at 1.1-km resolution, in terms of latitude and longitude. These values were used as the starting points for the ray casting during the reconstruction. The GMP provides the camera viewing zenith and azimuth values for each pixel at 1.1-km resolution. They were used to define the direction vectors of the cast rays.

The TC_CLOUD product provides cloud motion (speed) at 17.6-km resolution, in northward and eastward directions, and cloud-top heights, both with and without wind correction, at 1.1-km resolution [25]. All of the variables are reported at the reference time of the AN camera. The cloud motion variables were used to represent the wind field during the wind correction. The two variables for cloud-top height were combined such that, given the heights with wind correction, all their no-retrieval values were replaced with the corresponding heights without wind correction that have fewer no-retrieval values.

All of the products were interpolated as needed to match the pixel resolution of whichever of the two input cloud masks was used in the reconstruction. The details behind the interpolations can be found in Reference [31].

4.3. Custom Cloud Mask

By design, the RCCM attempts to be clear-sky conservative, which at 1.1-km resolution leads to significant overestimation of cloud fractions, particularly for typical trade-wind cumulus clouds [32,33]. As noted in Section 3.2, this could lead to significant overestimation of cloud volume, particularly for typical fair-weather cumulus sizes [34]. The RCCM alone, with its 1.1-km resolution, may be too coarse for many cloud scenes of interest. Hence, to demonstrate improved reconstructions with a higher resolution than the operational RCCM, we produced the custom cloud mask at 275-m resolution.

The custom cloud mask was derived from the RGB channels for a few oceanic scenes selected for reconstruction whose results are shown in Section 5. Two different tests were used in the pixel classifications and then combined: Red channel threshold and whiteness deviation threshold.

The red channel threshold test takes advantage of the fact that clouds in oceanic scenes almost always scatter more solar radiation than the background (except in occasions of sun glint), especially at the red channel. It puts a threshold on the red channel’s scaled radiances (obtained by right-shifting Radiance/RDQI by 2 bits) within the scene, similar to the method explored in References [32,34]. In other words, for every pixel with the red channel scaled radiance ${\mathit{R}}_{\mathit{r},\mathit{i}}$, where ‘$\mathit{r}$’ represents ‘red band’ and ‘$\mathit{i}$’ the pixel, it is classified as cloudy if ${\mathit{R}}_{\mathit{r},\mathit{i}}\ge {\mathit{T}}_{\mathit{R}\mathit{C}}$, where ${\mathit{T}}_{\mathit{R}\mathit{C}}$ is the threshold value.

The whiteness deviation test takes advantage of the fact that clouds are “whiter” in color than the oceanic background. It puts a threshold on the largest of the differences between each of RGB radiances (obtained by multiplying scale factor to scaled radiance) and their mean. In other words, for each pixel with the RGB radiances ${\mathit{L}}_{\mathit{r},\mathit{i}}$, ${\mathit{L}}_{\mathit{g},\mathit{i}}$, ${\mathit{L}}_{\mathit{b},\mathit{i}}$, and their average value ${\mathit{L}}_{\mathit{a}\mathit{v}\mathit{g},\mathit{i}}=\frac{{\mathit{L}}_{\mathit{r},\mathit{i}}+{\mathit{L}}_{\mathit{g},\mathit{i}}+{\mathit{L}}_{\mathit{b},\mathit{i}}}{\mathbf{3}}$, it is determined cloudy if $\frac{\mathbf{max}\left(\left|{\mathit{L}}_{\mathit{r},\mathit{i}}-{\mathit{L}}_{\mathit{a}\mathit{v}\mathit{g},\mathit{i}}\right|,\left|{\mathit{L}}_{\mathit{g},\mathit{i}}-{\mathit{L}}_{\mathit{a}\mathit{v}\mathit{g},\mathit{i}}\right|,\left|{\mathit{L}}_{\mathit{b},\mathit{i}}-{\mathit{L}}_{\mathit{a}\mathit{v}\mathit{g},\mathit{i}}\right|\right)}{{\mathit{L}}_{\mathit{a}\mathit{v}\mathit{g},\mathit{i}}}\le {\mathit{T}}_{\mathit{W}\mathit{D}}$, where $\mathbf{0}<{\mathit{T}}_{\mathit{W}\mathit{D}}<\mathbf{1}$ is the threshold value.

In both tests, the threshold values were chosen through human-eye trial-and-error tests with the RGB image. For different view angles, different values were set for the red channel threshold, while a single value was set for the whiteness deviation threshold.

Table 1. Threshold values used in the custom cloud masks for the selected scenes.

Scene		Orbit 19726 Block 86		Orbit 1181 Block 91
Threshold Test		Red Channel $\left({\mathit{T}}_{\mathit{R}\mathit{C}}\right)$	Whiteness Deviation $\left({\mathit{T}}_{\mathit{W}\mathit{D}}\right)$	Red Channel $\left({\mathit{T}}_{\mathit{R}\mathit{C}}\right)$	Whiteness Deviation $\left({\mathit{T}}_{\mathit{W}\mathit{D}}\right)$
View Angle	DF	1201	0.4	2069	0.3
	CF	1150	0.4	1702	0.3
	BF	1116	0.4	1488	0.3
	AF	1202	0.4	1554	0.3
	AN	1179	0.4	1443	0.3
	AA	1114	0.4	1588	0.3
	BA	1065	0.4	1663	0.3
	CA	1100	0.4	1835	0.3
	DA	1187	0.4	2068	0.3

lists the threshold values used in the scenes selected for reconstruction. In the end, the cloud masks derived from each of the two tests were combined such that, if a pixel is classified as cloudy in either test, then it is classified as cloudy in the final mask. Relying on the accuracy of the human in choosing the appropriate threshold values, no other validation was conducted for the custom cloud mask. More discussion on this point can be found in References [27,32,34,35]. Visual comparisons between the RGB image and the custom cloud mask are shown in Section 5.

4.4. Wind Correction

The reconstruction for a scene concerns an instant in time, so all its input cloud masks at different view angles need to be adjusted to the same time. However, the approximate 7-minute difference between the DF camera and the DA camera leads to possibly significant cloud motion and displacement by wind, as well as cloud morphing. Hence, there arises the need for wind correction to offset the cloud displacement and morphing between each view angle and the reference time of the reconstruction.

Since the cloud motion and height variables for each block are all reported at the reference time of AN (the nadir view angle), and because the maximum time difference between view angles is minimized when AN is the reference time as opposed to others (less than 4 minutes), all wind corrections and therefore reconstructions in this study were conducted with reference to the time of AN.

The wind correction poses a significant challenge. The reported cloud motion parameters only concern the cloud top boundary points observed for each pixel in the AN camera, so it is difficult to infer the overall wind field at different altitudes throughout the scene. That is, there is no direct information for the wind in other view angles, where the pixels of the same geolocations are looking through different portions of the atmosphere and different parts of the observed clouds. The reported wind at AN cannot simply be superimposed as the wind field on the other view-angle images that may be viewing parts of the cloud at a lower altitude.

With the above constraint in mind, a method named “cloud spreading” (“cloud thickening” in Reference [31]) has been devised. First, the greatest wind value is selected from among the reported wind values in the block. Then, regardless of the direction of the wind, its quantity (speed) is used to “spread” in all horizontal directions every cloudy pixel in the block. The level of spreading is determined by the time difference between the view angle and the AN, multiplied by the chosen maximum wind speed. Hence, more cloud spreading is applied to more oblique view angles. In other words, for each cloudy pixel ${\mathit{P}}_{{\mathit{i}}_{\mathbf{0}},{\mathit{j}}_{\mathbf{0}}}$, every other pixel in the block ${\mathit{P}}_{\mathit{i},\mathit{j}}$ is assigned ${\mathit{P}}_{\mathit{i},\mathit{j}}={\mathit{P}}_{{\mathit{i}}_{\mathbf{0}},{\mathit{j}}_{\mathbf{0}}}$, if its index numbers, $\mathit{i}$ and $\mathit{j}$, satisfy ${\left(\mathit{i}-{\mathit{i}}_{\mathbf{0}}\right)}^{\mathbf{2}}+{\left(\mathit{j}-{\mathit{j}}_{\mathbf{0}}\right)}^{\mathbf{2}}\le {\mathit{N}}^{\mathbf{2}}$, where $\mathit{N}=\mathit{r}\mathit{o}\mathit{u}\mathit{n}{\mathit{d}}_{\mathit{u}\mathit{p}}\left(\frac{\left|\mathit{V}\right|·\Delta \mathit{T}}{\mathit{R}}\right)+\mathbf{0.5}$ is the spreading addition number in pixel unit, $\mathit{V}$ the selected wind vector with the block’s highest speed, $\Delta \mathit{T}$ the time difference between an oblique camera and the AN, and $\mathit{R}$ the pixel resolution.

This method operates on the assumption that the greatest wind reported is the actual greatest wind in the entire scene, in any 3-D direction, and throughout the time interval between different view angles. This is a reasonable assumption, since in reality (1) wind is usually stronger in higher altitudes, and (2) wind field does not change much in the span of several minutes.

The basic principle behind this cloud spreading method is that all pixels in non-AN view angle images that are “potentially” cloudy at the time of AN are labeled as cloudy. This works because the nature of the reconstruction method “tolerates” the spreading of cloudy pixels. The reconstruction method assumes that a cloud should have been observed by all cameras and takes as cloudy only the cloudy voxels that overlap in all view angles. Hence, the cloudy pixels labeled by the spreading method can be simply considered as the candidates for the cloud’s location at the time of AN.

This wind correction method adds another source of uncertainty that tends to contribute a high bias to the reconstruction outcomes. The bias would depend on factors, such as the true 3-D wind field and cloud movements, and depending on the maximum wind speed applied, the gross overestimation of cloud volume is inevitable. Therefore, the reconstruction results both with and without the wind correction are shown and compared in this study, as a means of bounding potential errors caused by cloud motions over the time interval that MISR takes to observe the same scene.

4.5. Ray Casting Calculations

The final step in the reconstruction is the ray casting, where rays defined in terms of the MISR geometric and geographic parameters are cast through the pre-defined reconstruction domain, and voxels go through a series of calculations to be classified as cloudy or clear. The method of voxel boundary intersections is used.

Given the pre-determined reconstruction domain consisting of voxels on the WGS84 ellipsoid, the voxel boundaries are defined either as a plane of the form $\mathit{A}\mathit{x}+\mathit{B}\mathit{y}+\mathit{C}\mathit{z}=\mathit{D}$ or as an ellipsoid of the form $\frac{{\mathit{x}}^{\mathbf{2}}+{\mathit{y}}^{\mathbf{2}}}{{\left({\mathit{R}}_{\mathit{e}\mathit{q}}+\mathit{H}\right)}^{\mathbf{2}}}+\frac{{\mathit{z}}^{\mathbf{2}}}{{\left({\mathit{R}}_{\mathit{p}}+\mathit{H}\right)}^{\mathbf{2}}}=\mathbf{1}$, where $\mathit{A}$, $\mathit{B}$, $\mathit{C}$, and $\mathit{D}$ are constants, $\mathit{H}$ the voxel’s height from the surface, ${\mathit{R}}_{\mathit{e}\mathit{q}}$ the equatorial radius of WGS84, and ${\mathit{R}}_{\mathit{p}}$ the polar radius (see Section 4.1).

Next, rays are cast from the geolocations of cloudy pixels to the camera locations, and the intersection points between the cast rays and the voxel boundaries are calculated and recorded. Then, each intersection point is checked with respect to all other voxel boundaries in the reconstruction domain to see on which side of each of the voxel boundaries does it lie. This allows for pinpointing where the point lies within the reconstruction domain, in terms of its adjacent voxels. If any voxel in the entire domain has an intersection point in its boundary, then it is determined cloudy.

The number of cast rays per pixel in one dimension (SOM X or SOM Y axis) is determined as ${\mathit{N}}_{\mathit{C}\mathit{R}}=\mathit{r}\mathit{o}\mathit{u}\mathit{n}{\mathit{d}}_{\mathit{d}\mathit{o}\mathit{w}\mathit{n}}\left({\mathit{S}}_{\mathit{h}}·\frac{\mathbf{tan}\left(\mathbf{90}°-{\mathit{\theta }}_{\mathit{m}\mathit{a}\mathit{x}}\right)}{{\mathit{S}}_{\mathit{v}}}\right)+\mathbf{2}$, where ${\mathit{N}}_{\mathit{C}\mathit{R}}$ is the number of cast rays, ${\mathit{S}}_{\mathit{h}}$ the horizontal voxel resolution, ${\mathit{S}}_{\mathit{v}}$ the vertical voxel resolution, and ${\mathit{\theta }}_{\mathit{m}\mathit{a}\mathit{x}}$ the maximum viewing zenith in the block.

The addition of 2 is to ensure that there are at least two cast rays per pixel in one dimension and therefore four in total, one at each corner of the pixel. This ensures that no voxel in the reconstruction domain remains unpenetrated in between the travel paths of all cast rays.

Again, this process is repeated for all nine view angles, and only the voxels that have been determined cloudy by all nine angles are eventually classified as cloudy. For more discussion on this and other ray casting methods, see Reference [31].

Afterward, the cloud-top height (TC_CLOUD) is applied to reduce overestimation of cloud volume at the top of the reconstructed cloud (as mentioned in Section 2). This is done by casting rays from the AN cloud mask and calculating the intersection points only for the voxels at or below the cloud-top height for each pixel. Only the voxels that remain cloudy under this constraint are kept to produce the final reconstructed cloud volume.

5. Reconstruction from MISR Data: Results

The results presented here come from hand-picked scenes that show a good vertical development of clouds favorable for multi-angle volume reconstruction from MISR data. Our choices were guided by the assessment on the reconstruction method with respect to cloud geometry (Section 3.1). Validation of the results is difficult due to the lack of objective truth for real-world clouds, but we managed to choose a scene for which a cross-validation was possible (Section 5.1).

For each scene, a total of four types of results are shown as follows.

• Input: RCCM. Wind correction: On.
• Input: RCCM. Wind correction: Off.
• Input: Custom cloud mask. Wind correction: On.
• Input: Custom cloud mask. Wind correction: Off.

The reconstruction voxel resolution was set as 275 m in all cases, which was anticipated to produce small overestimation for the cloud sizes in the selected scenes, according to Figure 5. For the RCCM inputs, the level of high-confidence cloud was set to cloudy and all other classification levels set to clear. Hence, during the ray casting stage, rays were cast only from pixels classified as high-confidence cloud. In all of the figures showing the results, the reconstructed cloud volumes are plotted with rectangular cuboid voxels, but each voxel in reality represents the slightly bent cuboid described in Section 4.1. The top of the reconstruction domain was set to be 1 km above the block’s cloud-top height (Section 4.2).

Two scenes were chosen for reconstruction, one within Orbit 19726 Block 96 and the other within Orbit 1181 Block 91.

Table 2. Numerical results of the reconstructions. To simplify the volume calculations, each voxel was assumed to be a rectangular cuboid. The entries are in cubic meters.

Scene	Orbit 19726 Block 86		Orbit 1181 Block 91
Cloud Mask	RCCM	Custom Cloud Mask	RCCM	Custom Cloud Mask
With Wind Correction	6.5600 × 1012	4.6026 × 1012	2.3604 × 1011	1.4959 × 1011
Without Wind Correction	4.8882 × 1012	3.6264 × 1012	1.1133 × 1011	0.3157 × 1011

shows the numerical results of their reconstructed cloud volumes.

5.1. Orbit 19726 Block 86

The local region selected for reconstruction within MISR Block 86 of Orbit 19726 is shown in red boxes in Figure 7. It is over the Pacific Ocean far north of New Zealand near the Equator, observed on September 2, 2003, around 22:42:05 UTC. This scene was chosen for its cloud geometry favorable for reconstruction and to provide a cross-validation with the result from References [21,22].

Figure 7 also shows the RCCM and the custom cloud mask for the AN (nadir) view. The thresholds for the custom cloud mask in Table 1 were chosen specifically for the red box, but applied here over the entire block. The effects of sun glint on the right side of the image are often misclassified as cloudy, and the thinner clouds on the left side of the image are often classified as clear. This phenomenon is greater within the custom cloud mask. Except for the sun glint, the overall greater cloud cover in the RCCM is due to its coarser 1.1-km resolution as opposed to the custom cloud mask’s 275-m resolution.

Figure 8 shows for different image types the ray casting domain, which is slightly larger than the reconstruction region. The ray casting domain’s dimensions were determined by calculating the maximum possible horizontal travel distance of any cast ray from the largest viewing zenith of the entire block and the reconstruction domain’s height. This was to ensure that any pixel whose ray travels through the reconstruction domain is included in the ray casting domain. Figure 9 shows the reconstructed volumes with and without the wind correction (cloud spreading), from the RCCM and the custom cloud mask for the selected region. In all results, a tall and wide wall of cloud extending in the cross-track direction is seen.

The results are shown in Table 2 and visualized in Figure 9. Overall, the reconstructions from the RCCM exhibit larger volumes than those from the custom cloud mask, especially in the lower-altitude clouds that appear to be small cumuli in the RGB image. With clouds that are just a few hundred meters in size and observed by a coarser-resolution and clear-conservative mask like the RCCM, a large overestimation in cloud volume is expected as discussed in Section 3. Moreover, the greater cloudiness in the cloud masks with wind correction led to the greater reconstructed volumes than without wind correction. In all cases, the wall of clouds stretching from the upper left to the lower right corner of the reconstruction region is seen, as imaginable from our mental constructions from the nine RGB images.

The reconstruction from the custom cloud mask without wind correction for this scene is compared with the result from Reference [21] which also did not include wind correction, as shown in Figure 10. Reference [21] used a stereo-photogrammetric method to retrieve both cloud-top heights and sides with a precision of about 200-300 m. In terms of the general shape, the reconstruction agrees with their result, showing the wall of clouds that stretch vertically up to around 8.5 km. The reconstruction shows greater horizontal width between around 2 and 5 km in altitude, as [21]’s result takes a concave shape. When the observed cloud has concavity in the along-track direction, our reconstruction method would be unable to detect it and create a convex volume with greater overestimation.

5.2. Orbit 1181 Block 91

The second region selected for cloud volume reconstruction, shown in red boxes in Figure 11, is over the Atlantic Ocean south of West Africa, observed on March 8, 2000, around 11:33:19 UTC. This scene was chosen as it contains two vertically tall and horizontally narrow columns of clouds. As in the previous example, the thresholds (Table 1) were chosen specifically for the red box, but applied here over the entire block. Figure 12 shows the ray casting domain around the selected reconstruction region.

The results are shown in Table 2 and visualized in Figure 13. In all the reconstruction results, the semblance of two tall and tilted columns of clouds is seen. As in the case of Section 5.1, the RCCM contains more cloudy areas and therefore thicker reconstructed volumes. However, compared to Section 5.1’s case where the RCCM’s reconstructed volumes are up to about 1.4 times greater than the custom cloud mask’s, this section’s case shows up to 3.5 times greater. Besides the uncertainties in the custom cloud mask itself, such difference between the two cases could be attributed to the difference in the size of the observed clouds. The overestimation of cloudiness in cloud masks with coarser resolutions due to the resolution effect is relatively greater when the observed cloud is smaller in size.

The wind correction also exhibits greater reconstructed volumes than no wind correction. Compared to Section 5.1’s case where the wind correction results in reconstructed volumes about 1.3 times greater than without the wind correction, in this section’s case it does up to 4.7 times greater. Such difference between the two cases stems from the same difference as before, the size of the observed cloud. In general, if the observed cloud is smaller in size, it occupies fewer cloudy pixels in the satellite images. Hence, during the cloud spreading, under the condition of the same maximum wind speed applied, the ratio of the number of cloud-spread pixels to the original number of cloudy pixels in the cloud mask is greater. This leads to greater reconstructed cloud volumes. In this section’s case, the wind correction results in much greater volumes notwithstanding the fact that the maximum wind speed applied during the cloud spreading was smaller than Section 5.1’s case (10.5587 m/s as opposed to 14.1646 m/s).

6. Discussion and Conclusions

This paper explores 3-D cloud volume reconstruction from MISR using a ray casting method. First, simulations were done, with a handful of synthetically generated clouds and with respect to the observed cloud’s geometry and the reconstruction voxel resolution, to understand and estimate the errors and uncertainties associated with the reconstruction method in various cases. Under the assumption of perfect cloud detection, it has been demonstrated that overestimation of cloud volume in the reconstruction would be smaller when the observed cloud is more developed in the vertical and cross-track directions and less developed in the along-track direction. In addition, overestimation would be smaller when the size of the observed cloud is large and the reconstruction voxel resolution is small in comparison to the pixel resolution.

Then, reconstructions were done with the real MISR data using the RCCM with 1.1-km resolution and the custom cloud mask with 275-m resolution generated from MISR RGB radiance data. The Earth’s surface curvature was reflected in the reconstruction domain and wind correction was performed by implementing a method termed cloud spreading. Reconstruction results were shown for two hand-picked scenes, one within Orbit 19726 Block 86 and the other within Orbit 1181 Block 91. The former case showed some degree of resemblance with the result from Reference [21]. Several tendencies were observed in both scenes. First, the RCCM results tend to overestimate more than the custom cloud mask results. This is mainly because the greater ratio between the pixel resolution and the size of the observed cloud leads to greater overestimation of cloudiness in the cloud masks [32]. The greater precision (finer details) of the custom cloud mask results, and the error evaluation with the synthetic clouds in Section 3, point to a need for MISR to move toward a cloud mask product with 275-m resolution. Next, the results with the wind correction tend to overestimate more than the results without the wind correction, as anticipated in Section 4.4.

While we conclude that the reconstruction results are visually positive and agreeable to the human eye, and that 3-D cloud volume reconstruction is possible with MISR using the ray casting method and proper cloud masks, several limitations need to be noted:

• Uncertainties in the reconstructions are difficult to ascertain and quantify due to the lack of truth for comparison. Their dependence on the true underlying structures of observed clouds and on the clear-conservative nature of cloud masks are similar to discussions surrounding measurements of cloud fraction and cloud horizontal size distributions (e.g., References [27,35,36]). One advancement that could be made in this regard is a proper coordination between MISR and ground assets (e.g., scanning cloud radars) deemed to provide reliable estimates of cloud shapes and volumes. Another could be a better validation of the reconstruction method by using 3-D radiative transfer simulations imitating MISR measurements applied to 3-D cloud fields generated from cloud resolving models. This could better reflect the nuances of cloud masking and cloud motion on the cloud volume reconstruction (cf. Reference [37]).
• Cloud masks with higher spatial resolutions are simply better for cloud volume reconstruction. Provided that the issue of computation is well handled, data from instruments with higher resolution could significantly improve the reconstruction results. Future work could make use of ASTER onboard the same Terra satellite as MISR. ASTER looks at the same scene as MISR at the finer resolution of 15 m and with two view angles, and hence could help refine the reconstructed volumes and generate better results.
• The application was limited to cases of vertically well-developed cumulus clouds that had been manually selected. The automation and operationalization of the reconstruction program for all available MISR data remains a challenging task, which may be resolved when a technique for recognizing applicable scenes is developed, such as properly trained machine learning algorithms or artificial neural networks.

While the limitations are legitimate and require a significant future effort to be resolved, there are several practical implications that give reasons to move forward with 3-D cloud volume reconstruction:

• Provided that the reconstructions are verified with proper and enough validations, they could be used to study the 3-D cloud distribution of a select geographical location for a climatological period of time, thanks to MISR’s long-lasting operation and stable 10:30-am equator crossing time of its sun-synchronous orbit. Such studies cloud lead to understanding the clouds’ impact on the local geography’s radiation budget and further our understanding of the role of clouds in the Earth’s radiation budget.
• Once the reconstruction program is extensively used to generate and accumulate results for many different real-world cloud scenes, the collection of such data could be used to test the realism of the geometry of clouds produced in large-eddy simulated models and cloud resolving models.
• The reconstructions could be useful in providing the initial guess for future tomographic research involving MISR. Tomography is a potential tool in remote sensing retrievals of 3-D atmospheric profiles, including cloud microphysical properties, and as explained in Reference [38], the tomographic reconstructions could greatly benefit in accuracy and computation if the initial conditions, such as cloud boundaries, are specified.

Peering into the future of passive remote sensing satellite missions for cloud volume reconstruction, several recommendations are presented. First, multi-angle instruments on two or more platforms in orbital formation could be designed to view the same scene from multiple directions at the same time. This would eliminate the uncertainties associated with cloud motion that MISR faces. Second, adding a view angle more oblique than MISR’s DF and DA cameras (±70.5°) would reduce overestimation due to the prism building effect at cloud top and bottom (Section 2 and Section 3). While increasing the number of view angles would lead to greater accuracy, experiments not presented in this study suggested that there is little improvement in adding more angles between the current range of MISR view angles. Finally, the spatial resolution of future imagers is advised to be around 50 m. This resolution is recommended to improve cloud masking by reducing sub-pixel cloud effects and increase textural cloud discrimination (e.g., References [34,35]), to reduce overestimation of cloudy area to negligible amounts for real-world clouds (e.g., Reference [36]), to be able to estimate cloud-base height based on cloud shadow detection (e.g., Reference [39]), and to provide improved precision in cloud stereoscopic height retrieval to the order of 10 m using a similar algorithm as in Reference [25]. These recommendations would provide a highly favorable mission scenario for employing the ray casting method for future 3-D cloud volume reconstructions.