Stereophotogrammetry of clouds observed during T-REX

Stereophotogrammetric images collected during the Terrain-induced Rotor Experiment (T-REX), which took place in Owens Valley, California, in the spring of 2006, were used to track clouds and cloud fragments in space and time. We explore how photogrammetric data complements other instruments deployed during T-REX, and how it supports T-REX objectives to study the structure and dynamics of atmospheric lee waves and rotors. Algorithms for camera calibration, automatic feature matching, and 3D positioning of clouds were developed which enabled the study of cloud motion in highly turbulent mountain wave scenarios. The dynamic properties obtained with photogrammetric tools compare well with data collected by other T-REX instruments. In a mild mountain wave event, the whole life cycle of clouds moving through a lee wave crest was tracked in space and time showing upward and downward motion at the upstream and downstream side of the wave crest, respectively. During strong mountain wave events the steepening of the first lee wave as it developed into a hydraulic jump was tracked and quantified. Vertical cloud motion increased from ~2 m s -1 to 4 m s -1 and horizontal cloud motion decreased from 20 m s -1 to 16 m s -1 with the development of the hydraulic jump. Clouds at distinct vertical layers were tracked in other mountain wave events: moderate southerly flow was observed in the valley (~8 m s -1 ), westerly motion of the same magnitude at the Sierra Nevada mountain crest level, and westerlies with speeds of over 20 m s -1 at even higher altitudes.


Introduction
Stereophotogrammetry is the reconstruction of three dimensional (3D) coordinates of points on an object by matching image pixels on two photographs capturing the same scene from different viewing angles. The reconstructed points are then combined into socalled "models" of the 3D surface of the object that is visible to the cameras. The methodology features qualities which make it highly valuable for studying ephemeral meteorological phenomena such as clouds: The models have high temporal resolution, limited only by the maximum frame rate of the cameras, and spatial resolution limited by the pixel size of the cameras and the distance between the object and the cameras. These characteristics enable the tracking of the evolution of cloud components of very fine scales, even down to the single turret level, providing a wealth of detail unparalleled by any other methodology. Furthermore, since all that is needed are two high quality digital cameras and two GPS sensors, the data collection process is relatively simple and costefficient.
Stereophotogrammetry has many different applications in various fields such as archeology (e.g. Orengo, 2013), topographic mapping (e.g. Mason et al., 2000), engineering (e.g. Bascoul et al., 1993), architecture (e.g. Solveig et al., 1979), and even police investigations (e.g. Breen and Anderson, 1986), and a large variety of commercial and open source software exists to process stereo photographic images (see e.g. en.wikipedia.org/wiki/Comparison_of_photogrammetry_software). However, when applied to cloud images, the processing methodology faces unique challenges, which make the available software packages largely unsuitable as the accuracy of the resulting model strongly depends on adequate camera calibration (see Sec. 3.2 for details). Since errors in the position of the 3D model data increase exponentially with the distance from the cameras to the object of interest, camera calibration algorithms optimized for applications where the objects are much closer to the camera often do not accomplish the level of accuracy needed for cloud images. Another difficulty lies in the fact that the available algorithms for matching pixels from one image to those of the other image often rely on the texture (edges, ridges, corners, color contrast, etc.) of the captured objectcharacteristics rarely observed in cloud images. In recent years stereophotogrammetric software has been successfully developed by several researchers (e.g. Hu et al., 2009, Öktem et al., 2014, Romps and Öktem, 2018 but to our knowledge none of them are publicly available. Therefore, our first objective was to develop image processing software which is able to handle the unique challenges presented by cloud images. The images used in this study were collected during the Terrain-Induced Rotor Experiment (T-REX) which took place in the spring of 2006 in the lee of the Sierra Nevada mountain range in Owens Valley, California, USA (Fig. 1, Grubišić et al., 2008). The main objective of T-REX was to explore the structure and dynamics of atmospheric lee waves, rotors 1 , and their interaction with the boundary layer. A conceptual model of a lee wave with rotor development is shown in Fig. 2: Downslope air motion in the lee of a mountain range leads to the development of lee waves. A rotor forms under the first lee wave crest which results in counter flow near the surface. Cap and rotor clouds form over the mountain and rotor, respectively. The geography of Owens Valley is uniquely conducive to the development of lee waves. It is located just east of the north-south oriented Sierra Nevada which forms a natural barrier to the westerly winds coming from the Pacific Ocean. The eastern slopes of the mountain range facing the valley are extremely steep (the vertical distance between the valley floor and the mountain reaches >3000 m) which leads to frequent occurrence of lee waves: an ideal setup for T-REX.
Various different instruments were deployed during T-REX including aircraft, airborne and ground-based radars, lidars, sodars, wind profilers, surface sensors, soundings, and also two digital photo cameras. An extensive dataset was collected during 15 Intensive Observation Periods (IOPs) which lead to various new insights published in numerous studies (see special collection published by the American Meteorological Society journals.ametsoc.org/topic/t-rex). However, the stereophotogrammetric images have not yet been analyzed. Our second objective is therefore to explore how the photogrammetric data can complement the existing data sets and expand our knowledge of the dynamics in the highly turbulent rotor environment.

Data
Stereo images from five T-REX IOPs (6b, 8, 9, 11, and 12) were used in this study. The images were captured with two Canon 20D digital 8 Mpixel cameras. The locations of the cameras were measured with GPS sensors and varied slightly among different IOPs but were located towards the southern end of Owens Valley. The baseline between the cameras was 453 m for IOP 8 and 666 m for all other IOPs investigated in this study. The frame rate was 0.1 fps for IOPs 6 and 9 and 0.2 fps for all others, which correspond to temporal resolutions of 10 s and 5 s, respectively. To maximize alignment, the cameras were first pointed at each other and then turned 90° horizontally toward the target until they were parallel. Finally they were pointed up by ~30°.
To verify and compliment the stereo images we used several data sets from other instruments. Horizontal wind, temperature, and humidity data collected by dropsondes launched from the UK FAAM BAe-146 aircraft (UK Met Office (UKMO), 2007) were used to investigate the vertical composition of the atmosphere and wind speed at different altitudes. Horizontal and vertical wind speed data measured by the University of Wyoming King Air Research Aircraft (UWKA, University of Wyoming, 2011) provide details on the position and characteristics of lee waves. These airborne instruments were complemented by two 915 MHz boundary layer wind profiler radars (ISS2 (UCAR/NCAR -Earth Observing Laboratory, 2011a) and ISS-MAPR (UCAR/NCAR -Earth Observing Laboratory, 2011b)) located farther north in Owens Valley as well as data from the Doppler radar from the German Aerospace Center (DLR, German Aerospace Center (DLR) (2007)). NCEP North American Regional Reanalysis (NARR) data (National Centers for Environmental Prediction/National Weather Service/NOAA/U.S. Department of Commerce, 2005) show the large-scale synoptic background.
Theoretical background: The pinhole camera model To discuss the theoretical background of stereophotogrammetric image processing we first set up two coordinate systems: The camera coordinate system and the image coordinate system (Fig. 3). The camera coordinate system is a right-handed system with three dimensions and we choose it in such a way that the left camera is located at the origin at x=y=z=0 (Fig. 3a). The y axis points from the left to the right camera, i.e. along the baseline so that the right camera is located at x=0, y=d, z=0, where d is the length of the baseline. Since the cameras were aligned in such a way that their main sight lines are parallel to each other and perpendicular to the baseline (Sec. 2) we choose the z axis parallel to the main sight line of the left camera and the x axis pointing up and perpendicular to the y and z axes (Fig. 3a). The image coordinate system (or image plain) is two dimensional with the axes pointing along the width (axis a) and the height (axis b) of the image (Fig. 3b). Points (or pixels) in the image plane are denoted with A i =(a i ,b i ). Before the camera calibration (Sec. 3.2) the image plane is parallel to the x-y plane of the camera coordinate system.
When looking at just one camera and assuming a pinhole camera (i.e., ignoring effects like lens thickness), all rays of light (or sight lines) from a 3D object in the camera coordinate system pass through the center of a camera and are mapped onto a 2D film or digital device, i.e. the image plane (Fig. 3b). Following Hartley and Zisserman (2003) this projection of a 3D object in the camera coordinate system onto a 2D image in the image coordinate system therefore follows the laws of projective geometry where the 3D camera coordinates of an object X i =(x i ,y i ,z i ) relate to the 2D coordinates on the image A i =(a i ,b i ) via a projection matrix P (note that in actual application homogeneous coordinates are being used, see e.g. Öktem et al., 2014): (1) The projection matrix P is the product of a camera matrix C, a rotation matrix R, and a translation matrix T = . (2) The camera matrix C contains all parameters related to the camera itself where f is the focal length of the camera, k a,b is the pixel density (number of pixels per meter, i.e. 1/(camera resolution) in a and b directions in meters), sk is the skew of the image plane (angle between the principal axes a and b of the image plane), and pp a,b are the principal point image coordinates in a and b directions of the image. The rotation matrix R describes the rotation between the actual camera (i.e. the image plane) and the camera coordinate system, where α, β, and γ are the three Euler angles turning counterclockwise around the a, b, and a theoretical c axis which is perpendicular to a and b, respectively, when looking towards the origin (Fig. 3b). Due to the setup of our coordinate systems we first assume all angles to be zero. However, the angles will change during the camera calibration process described in Sec. 3.2. The rotation matrices around each individual axis are ].
The translation matrix takes the location of the cameras into account: where x cam , y cam , z cam are the coordinates of the cameras in the camera coordinate system. As already mentioned, in our case x camL =y camL =z camL =0 for the left camera and x camR =0, y camR =d, z camR =0 for the right camera (Fig. 3a). If we take both the left and the right camera into account, we are left with two equations When all parameters in the projection matrices are known and the image coordinates of an object (A L and A R ) have been identified in the images, the set of equations can be solved for X using e.g. single value decomposition (Hu et al., 2009): the 3D coordinates of the object in the camera coordinate system are calculated via basic triangulation.

Camera Calibration
As mentioned above, accurate 3D positioning of an object relies on the detailed knowledge of the projection matrix parameters. The importance of precise camera calibration will be explained in the following thought experiment: Imagine two cameras on a horizontal plane focusing on an object in the same plane which are positioned in such a way that the object and the cameras form a triangle with all three sides of roughly equal length. If one horizontal angle of one camera is assumed wrong by a tiny fraction the calculated position of the object will be either too far away or too close to the cameras. If the error in the angle is small, the difference in the distance will also be small. However, if the distance between the cameras and the object is very large compared to the baseline (i.e. the triangle is very pointed) a small error in one camera angle will have a large effect on the calculated position of the object. Since clouds are by nature far away from the cameras, the distance between the cameras and the clouds compared to the baseline can easily reach a factor of 10, making accurate camera calibration of utmost importance.

a) Lens distortion correction
When using modern cameras instead of a simple pinhole camera, the effect of lens distortion needs to be considered before the pinhole camera model can be applied. The two most common forms of lens distortion are barrel and pincushion distortion where the curvature of the photo lens makes lines that are straight in the real world appear to be bent outward or inward on the image, respectively. The distortions are approximately radially symmetric with the strongest effect at the edges of the image. The exact amount of distortion depends on the individual lens of the camera.
To correct for lens distortion the standard procedure is to calculate new image pixel coordinates based on some polynomial function whose coefficients depend on the properties of the lens (see e.g. Hugemann, 2010). The best way to get these coefficients is to take several images of an object with many rectangular straight lines (e.g. a checkerboard) from different angles and calculate the lens distortion coefficients directly. Unfortunately, at the time of image processing, the cameras used in T-REX were not available to us so the best option was to get the coefficients from the lensfun library (http://lensfun.sourceforge.net/). The lens distortion correction function from the imagemagick software package (http://www.imagemagick.org/) was then used to apply the correction to the images. Obviously, the coefficients from the lensfun database are not from the exact lenses that were used to take the images and because our images do not contain any straight lines we have limited knowledge on how well the correction is actually working. (For an estimate of the overall measurement error see Sec. 3.4) b) The projection matrix: translation and camera matrices As described in Sec. 3.1, the projection matrix consists of the camera, rotation, and translation matrices. The translation matrix contains the positions of the cameras (Eq. 8) which were measured by GPS whose error range was given as 0.5 m. No further correction of the camera position was applied as accurate camera calibration is much less sensitive to camera location than to camera angles. 2 The parameters of the camera matrix are skew (zero in our case), principal point coordinates, pixel density, and focal length (Eq. 3). The principal point coordinates were assumed to be the center of the image, i.e., at half image width (½ * 3504 pixels) and height (½ * 2336 pixels). The pixel density was calculated from the camera resolution (6.42 μm) and the focal length was 10 mm for all IOPs. Note that these parameter values likely have some errors, however, including them in the optimization process described below did not lead to better results and so we considered them fixed at the given values.
c) The projection matrix: rotation matrix Getting the camera angles in the rotation matrix right (Eq. 4) is the most essential but also most difficult part of the camera calibration process. To measure the deviation of the camera angles from the camera coordinate system we make use of a modified version of the "epipolar line constraint" (Hartley and Zisserman, 2003): In two dimensions, two sightlines from two cameras focusing on the same object always meet on the 2D plane they are located on, i.e., the triangulation of the object's position always has a solution no matter how erroneous the angles in the rotation matrices are. In three dimensions this is not necessarily the case. If the camera angles in the rotation matrices are wrong, the calculated sightlines will pass each other at a certain distance dd in 3D space (e.g. the magenta and dark blue lines in Fig. 4). Since we know that the sightlines do meet in reality (they are pointed at the same object) we can use this principle to calibrate the cameras by adjusting the camera angles in such a way, that the distance between the 2 Imagine two cameras with a baseline of d=666 m, and a cloud located straight ahead of the left camera at a distance of z=10 000 m. The angle between the baseline and the sightline of the right camera is then 86.2°. If we now assume an error of 0.1° in the right camera angle β, leading to an angle of atan(z/d)=86.3° the distance z will be calculated as z=tan(86.3°)*666 m=10 299 m leading to a cloud location that is wrong by 299 m. However, if we assume the correct angle of 86.2° but a baseline with a 0.5 m error, i.e. d=666.5, we calculate a distance z=tan(86.2°)*666.5 m=10 035 m leading to a cloud location that is wrong by only 35 m. sightlines is at a minimum (ideally zero): First we identify the image coordinates of an object in both images (pixel pairs a 1L , b 1L and a 1R , b 1R, Fig. 4). Then, for each camera separately, we assume initial camera angles (in our case all zero). Using equation 11 below, we calculate a vector v for each camera pointing from the camera position in the direction of the sightline between the camera and the object The shortest distance dd between the vectors v L and v R is where the line connecting the two sightlines is perpendicular to both sightlines, i.e., a line along the cross product between the two vectors n=v L ×v R . Since the location of the left camera is x camL =y camL =z camL =0 the shortest distance dd min between the sightlines is To calibrate all six camera angles, theoretically we need six pixel pairs to minimize the six resulting distances dd 1…6 (two examples are conceptually shown in Fig. 4). However, since the identification of the pixel pairs is subject to errors (see Sec. 3.5 on automatic feature matching) we use at least 70 pixel pairs for each IOP in order to find the optimal combination of camera angles while minimizing all dds. In practice we use the GenSA optimizer from the statistical software environment R (Xiang et al., 2013) to find the optimal camera angles in order to calculate the projection matrices (Eq. 2) for the triangulation (Eqs. 9, 10).

Georeferencing
After successful camera calibrations for each IOP, we are now able to calculate the 3D coordinates of an object in the camera coordinate system. As a next step, we need to transform each pixel from the camera coordinate system (x, y, z) to the geographic coordinate system (longitude, latitude, altitude above mean sea level). To transform a 3D coordinate system, three points with coordinates known in both coordinate systems are needed. The GPS coordinates of the cameras provide two such reference points and so theoretically we need just one more. We were able to identify between two and four mountain tops or ridges in the images for each IOP which could be associated with respective features in GoogleEarth which provides the geographic coordinates. Unfortunately, as the cameras were pointed up at the clouds, the reference points were all located at the lower edges of the images where the possible error in the lens distortion and its correction is likely the largest (Sec. 3.2. a).
Once the reference points with known locations in the geographic coordinate system have been identified on the images, their camera coordinates can be calculated by triangulation. The camera coordinate system is then transformed to the geographic coordinate system in such a way that the coordinates of the cameras are fixed and the distances between the reference point positions in the camera and geographic coordinate systems are minimized resulting in the georeferencing rotation matrix.

Error analysis
As mentioned in the previous subsections, most of the image processing steps are prone to inaccuracies and so it is necessary to analyze how much uncertainty is inherited in the results. To assess the error in the calculated positions of the objects, we used the same reference points as for georeferencing (Sec. 3.3). As the bulk of the error is expected to be the result of inaccurate camera angles and not erroneous georeferencing, we believe that using the same reference points is a valid approach.
The position error for each reference point E point was calculated as the difference between the reference points in camera and geographic coordinates (RR in Fig. 5) divided by the distance between the reference points in geographic coordinates and the left camera (x point in Fig. 5). As mentioned before (e.g. in footnote 2, Sec. 3.2. b), the distance between the object and the cameras is a nonlinear function of the camera angles. Hence, uncertainties in the camera angels will lead to a nonlinear dependence of the error on the distance. Therefore we fit a nonlinear power law to the reference point errors (E point ) vs distances (x point ) resulting in an overall error of E(x)=1.527*(x/km) 0.862 . The estimated error at the cloud-camera distances observed in the different IOPs (5 to 15 km) is therefore ~6-15%. We consider an error of this magnitude acceptable for our purpose, especially because each error is constant within its respective IOP and therefore possibly somewhat changes the quantity of our results but has little effect on the qualitative findings. The relatively small magnitude of the error is further confirmed by the excellent agreement between photogrammetrically calculated cloud parameters with data from other instruments as described in Sec. 4.

3.5.
Automatic feature matching To triangulate the position of an object (e.g., a cloud "feature") it is necessary to identify the exact pixel the object is captured on in both images. This feature matching is usually carried out in a two-step process: 1. Characteristic features are identified in one image. 2. These features are matched with the corresponding pixels in the second image.
As already mentioned, the identification of characteristic features relies heavily on the structure of the captured objects like edges, lines, color contrast, etc. However, cloud images are generally lacking the necessary structure and so commercial feature matching algorithms like e.g. Agisoft (www.agisoft.com) or the feature matching functions of OpenCV (opencv.org), an open source computer vision library (Bradski, 2000) fail already at step 1. Therefore, we developed new software which does not identify characteristic features but rather tests all pixels in a chosen cloud feature for possible matches as described below. This new software follows the approach of Öktem and Romps (2015) and also makes use of some OpenCV functions.
Since we are mostly interested in a small part of an image (e.g. a single cloud), we first apply a mask using a mean shift segmentation algorithm developed by Christoudias et al., 2002, and modified by Joe Wezorek (jwezorek.com) and the OpenCV mask utilities to mask out everything outside the region of interest. The masked left image with a marked starting pixel is then displayed on the screen and the user chooses an approximate match for the starting pixel in the right image. Two types of rectangles are now defined: a) One rectangle around the user defined point in the right image, i.e. a rectangle that likely contains the exact match (the search region). b) Two rectangles of equal size around the starting point in the left image and around the possible matching point in the right image (the correlation regions of the left and right image). For each point in the search region a cross correlation coefficient is calculated between the correlation regions of the left and right images. The point with the highest cross correlation coefficient is considered a match for the starting point in the left image if it exceeds a certain predefined threshold. The process is now repeated for the next pixel in the mask of the left image. However, since it can be assumed that the distance between the new left pixel and the new right pixel match is similar as the distance between the previous left pixel and matched right pixel, a point on the right image calculated from this previous distance can now serve as the approximate right pixel around which the search region is defined and no additional user input is required.
For large masks (i.e. large clouds) the assumption made above that the distance of the new pixel pair is similar to that of the old pixel pair is not necessarily true. Therefore, if the mask exceeds a certain size it is divided into sub-regions. The user now defines a starting point for each sub-region at the beginning of the feature matching procedure. This splitting into sub-regions also has the advantage that the software routine can be parallelized which greatly increases the performance speed. The resulting list of matched pixels is checked for duplicate matches, i.e., right image pixels that are matched to two different left image pixels. If such a duplicate match is found, the match with the higher correlation coefficient is chosen and the other one is discarded.
Depending on the structure of the cloud image and the selection of the correlation coefficient threshold, more or less matches are found which, after triangulation, result in higher or lower density 3D cloud models. An example of a cloud and the matched pixels is shown in Fig. 6a and b. The 3D positions of the points are calculated and plotted in Fig. 6c. The 3D cloud model in Fig. 6c is rotated and captured from different viewing angles to illustrate the reconstruction of the outer surface of the cloud with photogrammetry. With the georeferencing matrix calculated in Sec. 3.3 the 3D cloud model is transformed to longitude, latitude, altitude coordinates.
Because only two cameras were deployed in T-REX, the 3D cloud model is restricted to the "surface" area of the cloud that is visible to both cameras. The use of several stereo cameras pointed at the cloud from different viewing angles, as carried out by Romps and Öktem, 2018, can increase the modeled portion of the cloud.

Preparation and analysis of cloud data
In principle, the 3D cloud models generated in the previous section can be used for analysis of properties such as cloud growth, or cloud volume. However, in our application we are only interested in the motion of a cloud or cloud fragment as a whole. Therefore, we average the latitude, longitude, and altitude coordinates of all matched pixels in an image pair to get the coordinates of the center point of the cloud or cloud fragment as seen by the cameras. We repeat the whole process described in the previous sections and the calculation of the cloud center for every image in a temporal sequence so that we are left with one cloud location in geographic coordinates for each time step, which we will call the "cloud track".
To characterize the cloud motion for a specific T-REX IOP it is desirable to calculate average cloud tracks for all clouds observed in the same environment. In order to calculate these track averages we normalize each cloud track by the 4D coordinates at the first time step, i.e. we pretend that all cloud tracks start from the same location in space and time. Note that at this point the clouds still move in different directions from this common starting point. However, we can make use of two observations (Sec. 4): a) each individual cloud roughly moves along a straight line and b) clouds captured in the same environment generally move in parallel. Therefore, for each individual cloud track we fit a straight line to the data points, which is treated as one single horizontal coordinate along the cloud track. We then average over all fitted lines, assuming that clouds from the same environment move in parallel, which leaves us with only one common horizontal coordinate along the average cloud track. The horizontal orientation of this along track coordinate gives us the direction of the cloud motion in degrees. We use the meteorological convention where 0˚ describes motion from south to north, 90˚ motion from east to west, etc.
With the knowledge of both the geographic coordinates of the cloud center and the temporal resolution of the images, we are able to calculate the horizontal and vertical cloud speeds. The horizontal cloud speed and cloud motion direction are expected to be related to wind speed and direction at the cloud level. However, this relation is only approximate as the cloud speed is a combination of external wind speed and internal cloud processes that include e.g. turbulence, condensation and evaporation. We compare the cloud speeds to wind speeds measured with other T-REX instruments in the following section.

Results and discussion
4.1.
IOP 9: Life cycle of clouds in a lee wave crest IOP 9 took place on April 2, 2006, and, according to the T-REX Science Director Summary, "featured a relatively weak mountain wave event". Probably because of its moderate nature, to our knowledge, data from IOP 9 have not been analyzed and published in the literature. However, it is exactly the relatively straightforward nature of the event that makes it an ideal candidate to test our methodology. The synoptic conditions (Fig. 7a) show west-southwesterly winds over California at the leading edge of a trough approaching from the west. The winds over Owens Valley are perpendicular to the Sierra Nevada mountain range: an ideal set up for lee wave activity. We analyze data that were collected in the morning hours (8:30 to 10:30 Local Time (LT)). During that time frame, there were generally clear sky conditions (Fig. 8) but small cloud fragments formed in the lee of the Sierra Nevada Mountains, presumably at the location of the first wave crest. In general, clouds formed on the western side of the wave crest and grew to larger sizes as they moved up into the crest where they remained stationary and often merged with other clouds. However, some of the clouds did not grow to significant scales and were carried with the west-southwesterly winds through the wave crest and down until their dissipation. By focusing on these small scale clouds we were able to track their location throughout their life cycle. Twenty clouds were tracked and processed as described in Sec. 3.6 and the results are shown in Fig. 9. Generally, the clouds' life cycle was completed in ~3 minutes as they moved ~2.5 km (Fig. 9a). The average horizontal cloud speed (Fig. 9c) was around 15 m s -1 (Table 1) moving from west-southwest to east-northeast (Table 1). The clouds formed at an altitude of ~4500 m and moved up into the first wave crest with positive vertical speed and then down until their dissipation ( Fig. 9b and d). The horizontal tracks of the individual clouds are shown in Fig. 10a and b (yellow lines).
Vertical wind data collected from the WY King Air aircraft during the same time frame are shown in Fig. 10c as the plane moves along the track shown in Fig. 10a and b. The upward (red) and downward (blue) air motion in the first and second wave crests are clearly visible. The flight leg visible in Fig. 10c just above 4500 m (i.e. the level of cloud formation) is also shown on a horizontal map in Fig. 10a. It is striking how well the location of the first wave crest (i.e., the transition from upward to downward motion) recorded during this flight leg agrees with the location of the clouds when viewed relative to the terrain. Horizontal wind speed and direction from data from the same flight leg are shown in Fig. 10b. The wind speed is with ~15 m s -1 the same as the cloud speed (Table 1). The wind direction vectors shown along the flight track are parallel to the cloud tracks confirming the measured cloud direction. In addition, Fig. 10b and d show the location and horizontal wind speed data from two dropsondes launched from the Bae146 aircraft. The location of the sounding at cloud level is indicated in red in Fig. 10b. The track of the sonde again parallels the cloud tracks and the wind speed measured by the sonde at ~4500 m is also ~15 m s -1 (Fig. 10d).
The excellent qualitative and quantitative agreement between the data sets from different instruments and the cloud data in the relatively simple setup of IOP 9 confirms the quality of the photogrammetric data processing methodology and the value of the derived data.

4.2.
IOPs 6 and 11: Dynamics at the leading edge of rotor clouds IOP 6 took place on March 25, 2006, and was characterized by the approach and passage of a depression and associated front from west to east over Owens Valley (see e.g. Sheridan and Vosper, 2012, Fig. 7b). Winds were strong from a west-southwest direction and increased with the frontal passage to become the strongest wind storm observed during the whole field campaign. IOP 6 features one of the most pronounced rotor developments observed during T-REX and has therefore been studied extensively in all its different aspects (e.g. Grubišić et al., 2008, De Wekker and Mayor, 2009, Jiang and Doyle, 2009, Hill et al., 2010, Cohn et al., 2011, Sheridan and Vosper, 2012, Kühnlein et al., 2013, Strauss et al., 2016 so we aim to integrate the photogrammetric analysis with the existing literature. Synoptic conditions before the passage of a frontal system are known to be conducive to the development of rotors (Holmboe and Klieforth, 1957) and indeed, a well developed rotor cloud can be observed in the stereo images captured in the morning hours before the frontal passage (Fig. 11a). IOP 11, which took place on April 9, features another case of strong rotor development which attracted a lot of interest from the scientific community (e.g. Drechsel et al., 2009, Armi andMayr, 2011). The large-scale synoptic conditions were somewhat similar to those observed during IOP 6, again featuring southwesterly winds at the leading edge of a trough approaching from the west with an embedded frontal system (Armi and Mayr, 2011, Fig. 7c). Armi and Mayr (2011) give a very detailed description of the temporal evolution of the conditions in Owens Valley before, during and after the time frame of the frontal passage in the afternoon hours. They found that during the earlier hours the temperature and wind evolution resembles what one would expect for a thermally induced valley wind system. However, starting at around 15:30 LT large scale advected cold air starts to overflow the Sierra Nevada Mountains and westerly flow penetrates deep into the valley leading to the development of a hydraulic jump. Both phases i.e. conditions without and with a hydraulic jump (phase A and B, respectively) can be clearly distinguished in the stereo images ( Fig. 11b and c). Both phases feature prominent rotor clouds, however, the cloud in phase A appears thinner and more fractured than in phase B. Fig. 11: Images captured during IOP 6 (a) and IOP 11, phase A (b) and phase B (c).
Rotor clouds in IOP 6 and IOP 11 phase A and B were largely stationary which made them undesirable candidates for tracking. However, small cloud fragments developed just upstream of the main rotor clouds and are subsequently advected to and ultimately merged with the main clouds. We focus our attention on these upstream cloud elements to gain an understanding of the dynamics at the leading edge of a rotor system in different phases of rotor development. We were able to track 37 cloud fragments in IOP 6 between 8:30 and 11:00 LT, and 20 each in phases [14][15]:20 LT) of IOP 11 (Fig. 12). Interestingly, we were able to track the clouds in IOP 6 for almost twice as long as clouds in IOP 11 (60 vs. 35 seconds, respectively) implying a formation farther upstream of the edge of the main cloud in IOP 6 (horizontal wind speeds were similar, see below). Fortunately, the shutter speed was doubled for IOP 11 leading to a similar number of total data points. Fig. 12: As in Fig. 9 but for clouds tracked in IOP 6 (red), and IOP 11 phases A (blue) and phase B (red).
Both horizontal and vertical cloud speeds are higher than in IOP 9 (Table 1). While the higher horizontal wind speed can probably be attributed to synoptic forcing (Fig. 7), the higher vertical wind speed likely reflects the fact that in IOP 9 the clouds are tracked at the very top of the wave crest whereas in IOPs 6 and 11 the clouds are tracked farther upstream at the ascending part of the wave. Also, the wave was likely less steep in IOP 9. As in IOP 9, the vertical wind speed decreases as the clouds approach the main rotor cloud at the top of the wave (cf. Fig. 9 and Fig. 12). Compared to clouds tracked in phase A, clouds in phase B of IOP 11 have lower horizontal (~20 m s -1 vs. 16 m s -1 , respectively, Fig. 12a and c) but higher vertical wind speeds (~2 m s -1 vs. 3.5 m s -1 , Fig. 12b and d) resulting in a steeper and faster ascent of the cloud fragments into the main cloud during the presence of a hydraulic jump. This result is consistent with the widely accepted theoretical model of the dynamics of a hydraulic jump as presented e.g. in the modeling study of Hertenstein and Kuettner (2005). The streamlines shown in their Fig. 8 are significantly steeper at the upstream side of the wave in the case of a hydraulic jump than those shown for a case without a hydraulic jump. Steeper streamlines imply a steeper ascent of potential cloud fragments with higher vertical and lower horizontal wind speeds, as observed in the photogrammetric data. The results for IOP 6 are in between those for the two phases of IOP 11 but leaning more towards the hydraulic jump phase (Fig. 12). The possible presence of a hydraulic jump in IOP 6 underlines the analysis of Kühnlein et al. (2013) who found "a rapid deceleration of the westerly downslope flow ending in a nearly vertical convergence line which resembles an internal hydraulic jump" during the time period of IOP 6 analyzed in the photogrammetric data.
The positions of the waves in IOP 6 have been extensively studied. Sheridan and Vosper (2012) and Strauss et al. (2016) show interpolated vertical velocity data from all three airplanes deployed during IOP 6. Both studies show at least two consecutive waves extending vertically through the entire troposphere, the first of which is located directly over Owens Valley. Cohn et al. (2011) analyzed the wave structure by applying a fit algorithm to data from three wind profilers which were aligned perpendicular to the Sierra Nevada on the valley floor (yellow diamonds in Fig. 13a). They place the first wave crest at the location of the middle wind profiler at the time when the stereo images were collected. This proposed wave location agrees very well with vertical wind speeds measured by the King Air aircraft at ~3700 m (Fig. 13a) which were collected just below the level of the tracked clouds (~4000 m, Table 1). When viewed relative to the terrain, the cloud tracks further confirm the wave location. They are located in the updraft region and merge with the main cloud upstream of the main rotor cloud which is likely located at the wave crest. As in IOP 9, the horizontal wind speed (~18 m s -1 ) and direction (243°) from the King Air flight leg (Fig. 13b) and a dropsonde launched from Bae146 ( Fig. 13e and f) agree well with the wind speed and direction calculated from the tracked clouds (Table 1, Fig. 12a and c). A plan position indicator (PPI) lidar scan (Fig. 13d) from the DLR Doppler Lidar (German Aerospace Center (DLR), 2011) which was located at the location of the green diamond in Fig. 13a shows the base of the main cloud at ~4000 m, the same altitude as the tracked clouds (Table 1), confirming not only the horizontal but also the vertical position of the clouds. The relative humidity data recorded by the dropsonde also shows a moist layer at the same level (Fig. 13f).
The altitude where the tracked clouds are initiated is significantly higher in IOP 11 (4885 m) than in IOP 6 (3972 m, Table 1). Neither aircraft nor dropsonde data are available for the level of the IOP 11 clouds. However, vertical wind speed data from a WY King Air flight leg at the higher altitude of ~6400 m show a clear wave pattern (Fig.  14) which possibly extends down to the cloud level of ~4800 m. A vertical reflectivity cross section of King Air cloud radar data collected at the same flight leg is presented by Armi and Mayr (2011). It places the cloud base of the western edge of the jump cloud at ~4800 m, which agrees with our findings of 4885 m (Table 1). It is interesting to note that the western edge of the rotor cloud, represented by the eastern end of the cloud tracks in Fig. 14, moves west during the transition from the time frame without (yellow lines) to the time frame with a hydraulic jump (red lines). This shift in location is likely related to the steepening of the wave and the resulting lower horizontal wind speeds mentioned above.

IOPs 8 and 12: Vertical wind shear
Because the main T-REX objective was to investigate mountain waves, most IOPs took place in conditions with a trough just off the coast of North America, leading to mid to upper level southwesterly winds perpendicular to the Sierra Nevada mountain range (Fig. 7). These conditions are also conducive to strong low level southerly winds in Owens Valley (Zhong et al., 2008) as surface lows often form north of Owens Valley ahead of the upper level trough axis resulting in a pressure gradient within the valley. These low level southerly winds have been reported in numerous T-REX IOPs (see e.g. the overview by Strauss et al., 2016) and their evolution has been studied numerically by Schmidli et al. (2011). IOPs 8 and 12 took place on April 1 and 11, respectively, under these common T-REX conditions (Fig. 7d and e) but with somewhat weaker southwesterly winds than e.g. in IOPs 6 ( Fig. 7b) described in the previous section. Jiang and Doyle (2009) characterize IOPs 8 and 12 as "moist IOPs" based on upstream soundings and subsequently state that "moisture in general tends to weaken mountain waves" even though this weakening effect seems to be less pronounced in IOP 8. Both, the weaker southwesterly winds and the weakening through moisture probably lead to the classification of IOPs 8 and 12 as "weak wave events" in the respective Science Director Summaries and no observed wave clouds in the stereo images (Fig. 15). However, they both show a unique feature: clouds formed not only over the Sierra Nevada but also in the low level southerly flow which gives us the opportunity to investigate the vertical wind shear between the above-mountain-crest and the in-valley levels. In IOP 8, which took place on March 31, the cap cloud over the Sierra Nevada mountains extends east over Owens Valley (marked as "Cross valley clouds" in Fig.  15a). It is largely stationary but small cloud fragments detach from the eastern edge of the cap cloud and between 15:00 and 16:10 LT we were able to track 30 of them in the stereo images until their dissipation (red dots in Fig. 16). As expected, they move from west to east at a horizontal cloud speed of ~9 m s -1 (Fig. 16c). The vertical cloud speeds are negative as the clouds move with the downslope winds until their dissipation (Fig. 16d). Clouds carried by the southerly flow in the valley enter the images at the top edge (marked as "Along valley clouds" in Fig. 15a) and 22 of them could be tracked during the same time frame as the cap cloud fragments until they move in front of the cap cloud from the cameras' perspective and are no longer distinguishable in the images. Their direction is from south to north (Table 1) at a constant horizontal speed of ~8 m s -1 (Fig.  17c). The vertical speed of these along-valley clouds is around zero as they move on a constant level of ~2700 m ( Fig. 17b and d). The setup is similar in IOP 12, but the cap cloud over the Sierra Nevada does not extend as far east (cf. Fig. 15a and b) and is confined to the areas above higher terrain (Fig. 18b). Unfortunately, the more western location of the cap cloud puts it at an unfavorable angle for tracking since the cloud movement is mostly directly towards the cameras, which makes positioning of the clouds difficult. In addition, few cloud fragments detach from the cap cloud in IOP 12 and the ones that do are not necessarily located at the eastern edge of the cap cloud but are often found at the top or the bottom of the cap cloud. Often cloud fragments even form and then dissipate in the vicinity of the cap cloud making it almost impossible to find a statistically significant sample of one type of cloud fragment which can be tracked and analyzed. Nonetheless, we tracked 21 cloud fragments between 9:00 and 10:00 LT which are associated with the cap cloud but their erratic motion is clearly reflected in the cloud tracks shown in Fig. 18b (green lines) and the statistics analysis presented in Fig. 16. They show very little horizontal movement with a cloud speed of less than 5 m s -1 (Fig 13c) and no vertical cloud speed is detected (Fig. 16b). The wind direction, even though generally from the southwest, has a large standard deviation of ±30° (Table 1). Fig. 17: As in Fig. 9 but for along valley clouds tracked in IOP 8 (red) and IOP 12 (black).
However, in IOP 12 there is an additional layer of clouds captured by the cameras during the above mentioned time frame which moves from west to east through the images (marked as "High cross valley clouds" in Fig. 15b), often appearing at the left edge of the image and either dissipating within the image or moving out of the image at the right edge. These clouds are advected in by higher-level winds above the cap cloud at ~4000 m ( Table 1) and 30 of them could be tracked and analyzed with photogrammetric tools. They move at high horizontal wind speeds of just over 20 m s -1 (Fig. 16c) but interestingly, they do not move down as the cap cloud fragments in IOP 8 but slightly upwards with a vertical cloud speed of ~1-2 m s -1 (Fig. 16b). Characteristics of the 21 invalley clouds tracked in IOP 12 are almost identical to the ones from IOP 8 (Fig. 17). They move at almost the same altitude, direction, and horizontal cloud speed (Table 1). The vertical cloud speed of the in-valley clouds in IOP 12 is also around 0 m s -1 but, contrary to the clouds from IOP 8, which seem to move at a slightly negative vertical cloud speed, the clouds from IOP 12 seem to move slightly up ( Fig. 17b and d). Because of the clouds' motion parallel to the main camera axes this statistically insignificant down or upward motion may be a result of small errors in the camera angle calibrations. The yellow and blue diamonds in (b) mark the location of the ISS-MAPR and ISS2 wind profilers, respectively. The dots in the dropsonde track in (a) mark the location where the dropsonde was at the altitude of the cross valley (red) and along valley (yellow) clouds. Fig. 18a shows the path of two dropsondes launched from the BAe146 aircraft during IOP 8. The sounding locations at the times when the sonde was at cloud level are indicated in yellow and red for the along-and cross-valley clouds, respectively. At higher altitudes the sounding moves from a southwesterly direction but the direction changes to westerly just at the altitude of the cap cloud, agreeing well with the directions of the tracked cross valley clouds (red lines, Fig. 18a). The horizontal wind speed of just under 10 m s -1 measured by the soundings at cloud altitude (~3500 m, Fig. 18c) compares well with the 9 m s -1 cloud speed calculated from the cloud tracks (Fig. 16c, Table 1). At low altitudes, the paths of the soundings turn north, paralleling the tracks of the along-valley clouds with good agreement in cloud speed and direction between the sounding data and the cloud track statistics (Fig. 17, Fig. 18a (yellow lines), Table 1). The relative humidity measured by the dropsondes (Fig. 18d) indicates a moist layer at the cap cloud level (~3500 m) which is especially pronounced in the western sounding. On close inspection, the eastern sounding shows an additional thin moist layer at ~2200 m which may be associated with the along-valley clouds. The lower altitude of this moist layer compared to the altitude of the along-valley clouds (~2700 m) could possibly be attributed to the fact that the sounding is located farther east than the clouds when viewed relative to the terrain, which would be consistent with a west to east slope of the moist layer within the down slope conditions indicated by the cloud tracks of the cap clouds (Fig. 16b).
No dropsondes were launched in IOP 12, however, data from two partial flight legs from the WY King Air aircraft collected at the level of the along-valley clouds (Table 1) at altitudes between 2000 m and 3000 m clearly show the southerly flow (pink arrows in Fig. 18b) indicated by the cloud tracks (red lines in Fig. 18b). The horizontal wind speeds from the western parts of the flight legs, which are located at the same distance to the terrain as the along-valley clouds, agree well with the horizontal cloud speeds of ~8 m s -1 (Fig. 17c). Data from two wind profilers located north of the tracked clouds near the flight path (blue and yellow diamonds in Fig. 18b) complement the data from the flight legs and the cloud tracks (Fig. 19). Both the western ISS2 and the eastern ISS-MAPR profiler show a clear shift in wind direction at ~3200 m ( Fig. 19b and e). Wind direction measured by the profilers at the level of the upper level (~4000 m) and along-valley clouds (~2600 m) are ~230-240° and ~160-190°, respectively, and agree well with the respective directions of 239° and 157° calculated from the cloud tracks (Table 1). The upper cloud level (~4000 m) wind speed measured by the western profiler (Fig. 19a) is with ~15 m s -1 somewhat lower than the respective cloud speed of ~20 m s -1 (Fig. 16c) but the wind speed at the level of the along-valley clouds agrees well with the ~8 m s -1 speed calculated from the cloud tracks. The lower level horizontal wind speed of ~15 m s -1 measured at the eastern profiler (Fig. 19d) is significantly higher than the cloud speed of ~8 m s -1 confirming increasing wind speeds towards the middle of Owens Valley as measured by the WY King Air (Fig. 18b). The vertical wind speed measured by the wind profilers is small (less than 2.5 m s -1 , Fig. 19c and f) and changes frequently from positive to negative indicating no organized wave pattern. In this variable environment, the upward motion of the higher-level clouds (Fig. 16b) can neither be validated nor invalidated but does not seem unrealistic. Note that the cloud altitude of ~2800 m of the cloud fragments associated with the cap cloud in IOP 12 (Table 1) are just below the transition zone between westerly and southerly winds ( Fig. 19b and d) which might help explain their unorganized behavior.

Conclusions
Stereo images collected during the T-REX field campaign were used to a) develop algorithms for camera calibration, automatic feature matching, and 3D positioning of clouds and b) study the motion of clouds in and around rotor systems. Stereo photogrammetric analysis of cloud images provides several unique challenges: The unfavorable geometry of a short baseline compared to the distance between the cameras and the clouds makes accurate camera calibration extremely important but the lack of reference points makes reaching this goal difficult. In addition, the blurry nature of clouds and resulting shortage of distinct features makes camera calibration even harder and traditional automatic feature matching strategies unsuitable. Using a combination of well established computer vision techniques and newly developed methods we created a set of software tools which is able to overcome these challenges and successfully reprojects cloud locations from 2D images into 3D space. Clouds and cloud fragments were tracked in space and time in 5 different T-REX IOPs and vertical and horizontal cloud speed as well as direction of motion were calculated. The results were compared with horizontal and vertical wind speeds and directions measured by multiple other instruments and the very good agreement between the results clearly demonstrates the potential of stereo photogrammetry in the atmospheric sciences.
Data from different IOPs were used to investigate different phenomena. In IOP 9 we track small clouds which form as air ascends to the top of a lee wave crest and dissipate on the descent at the downstream side of the wave crest. The small scale and ephemeral nature of the clouds provides the opportunity to study their whole life cycle. Analysis of the cloud tracks quantifies the ascent and descent of the clouds as well as their horizontal motion, and positions the first wave crest downstream of the Sierra Nevada at the western side of Owens Valley. The cloud parameters such as location, and horizontal and vertical speed compare well with data from the WY King Air and dropsonde data, underlining the quality of the photogrammetric analysis.
In contrast to the moderate wave activity observed in IOP 9, IOPs 6 and 11 show development of a mature rotor system with a pronounced rotor cloud forming at the top of the first wave crest. We tracked small cloud fragments which formed upstream of the main rotor cloud and were subsequently advected and ultimately merged with the rotor cloud. Analysis of the tracked cloud paths shows a clear distinction between the initial phase of IOP 11 and the following mature phase, which features the formation of a hydraulic jump as described in the literature. Increasing vertical but decreasing horizontal wind speed from the initial to the hydraulic jump phase points to a steepening of the wave as predicted by theory and modeling studies of hydraulic jumps. Cloud tracks from IOP 6 show similar characteristics as the clouds tracked in the hydraulic jump phase of IOP 11. As in IOP 9, horizontal and vertical cloud speeds and positions compare well with and complement the respective characteristics derived from other data sources.
IOPs 8 and 12 occur under less extreme conditions than IOPs 6 and 11 and little wave activity is detected. The stereo images show no rotor cloud but rather a cap cloud over the Sierra Nevada mountain range. The distinguishing feature of these two IOPs is the formation of clouds in low level southerly flow within Owens Valley. Clouds embedded in the westerlies above the ridge level of the Sierra Nevada as well as clouds carried by the low level southerlies were tracked and the resulting cloud characteristics again compare well with data from other instruments. The occurrence of the low level clouds points to the existence of a horizontal moist layer which is not well represented in the sounding data.
As mentioned above, the goal of our study was twofold: a) algorithm development and verification of the results and b) exploration of the potential of stereo photogrammetry in the investigation of phenomena in highly turbulent environments. In the current study we limit the data to a statistical approach, averaging not only over each individual cloud snapshot (by calculating the center of gravity) but also spatially over several clouds. The good agreement between these averages and other data sources shows the accuracy of the methodology and its great promise for future applications. Looking at averages already indicates new insights into the phenomena associated with rotor clouds. However, now that the methodology has been established, all available information without averaging could be utilized. By maintaining the information of each individual tracked cloud, properties such as horizontal and vertical cloud direction and speed (possibly as a proxy for wind speed) could be attained in great spatial and temporal detail. In addition, the 3D cloud models could be used to calculate physical properties such as cloud growth, volume, and shape, among others. Given the relatively simple and inexpensive method of data collection it does not seem too far fetched to imagine deployment of several stereo cameras to either provide entire 3D cloud models (as opposed to the partial cloud models used in this study which only represent the part of the cloud that is visible to the cameras) or larger spatial coverage and data collection by utilizing several camera pairs.