Contrail study with ground-based cameras

Photogrammetric methods and analysis results for contrails observed with wide-angle cameras are described. Four cameras of two different types (view angle <90 or whole-sky imager) at the ground at various positions are used to track contrails and to derive their altitude, width, and horizontal speed. Camera models for both types are described to derive the observation angles for given image coordinates and their inverse. The models are calibrated with sightings of the Sun, the Moon and a few bright stars. The methods are applied and tested in a case study. Four persistent contrails crossing each other, together with a short-lived one, are observed with the cameras. Vertical and horizontal positions of the contrails are determined from the camera images to an accuracy of better than 230 m and horizontal speed to 0.2 m s−1. With this information, the aircraft causing the contrails are identified by comparison to traffic waypoint data. The observations are compared with synthetic camera pictures of contrails simulated with the contrail prediction model CoCiP, a Lagrangian model using air traffic movement data and numerical weather prediction (NWP) data as input. The results provide tests for the NWP and contrail models. The cameras show spreading and thickening contrails, suggesting ice-supersaturation in the ambient air. The ice-supersaturated layer is found thicker and more humid in this case than predicted by the NWP model used. The simulated and observed contrail positions agree up to differences caused by uncertain wind data. The contrail widths, which depend on wake vortex spreading, ambient shear and turbulence, were partly wider than simulated.


Introduction
Contrails are linear clouds often visible to ground observers behind cruising aircraft.The conditions under which contrails form (at temperatures below the Schmidt-Appleman criterion) and persist (at ambient humidity exceeding ice saturation) are well known (Schumann, 1996).The dynamics of young contrails depend on aircraft emissions and wake properties and the pattern of contrails changes over their lifetime (Lewellen and Lewellen, 1996;Sassen, 1997;Mannstein et al., 1999;Jeßberger et al., 2013).Though contrails have been investigated for some time, still little is known about the full life cycle of individual contrails (Mannstein and Schumann, 2005;Heymsfield et al., 2010;Unterstrasser and Sölch, 2010;Graf et al., 2012;Minnis et al., 2013;Schumann and Graf, 2013).
Ground-based contrail observations may help to understand contrail dynamics and ice formation (Freudenthaler et al., 1995;Immler et al., 2008).The observable contrail cover can be predicted with contrail simulation models (Stuefer et al., 2005;Duda et al., 2009;Schumann, 2012).Such models require numerical weather prediction (NWP) data and traffic data as input and observations for validation.Here, we use NWP data as available from the Integrated Forecast System (IFS) model of the European Centre for Medium-Range Forecasts (ECMWF).Information on air traffic can, as of a few years ago, be received online from so-called flight radar data in the internet, including aircraft positions transmitted by Automatic Dependent Surveillance-Broadcast (ADSB), at least from the majority of aircraft which have such equipment (Jackson et al., 2005;de Leege et al., 2012).

U. Schumann et al.: Contrail camera observations
Wide-angle digital cameras have been used before to observe clouds (e.g., Seiz et al., 2007) and contrails (Sassen, 1997;Feister and Shields, 2005;Stuefer et al., 2005;Atlas and Wang, 2010;Feister et al., 2010;Mannstein et al., 2010;Shields et al., 2013).Whole-sky imagers using a fisheye lens image the full sky down, or nearly down, to the horizon (Shields et al., 2013).More narrow wide-angle cameras cover only part of the sky but with higher resolution.Besides color cameras, multispectral cameras recording in several spectral wavebands are also available (Feister and Shields, 2005;Seiz et al., 2007;Shields et al., 2013).Camera observations often reveal interesting cloud properties, but a single camera is insufficient to determine the distance of an observable object (LeMone et al., 2013).A video scene from a single camera allows determining the angular but not the linear cloud speed.Horizontal contrail positions can be estimated if the contrail altitude is known from other sources (Atlas and Wang, 2010).A network of cameras has been used for observations of upper atmosphere clouds (Baumgarten et al., 2009) and other objects (Shields et al., 2013).
In regions with dense air traffic, the sky is often full of contrails, and the assignment of individual observed contrails to specific aircraft requires accurate contrail altitudes besides traffic information.The analysis of aged contrails requires the trajectory analysis from aircraft flight routes to contrail positions.
Here, we report on a case study where we observed contrails with four wide-angle cameras, placed several kilometers apart and oriented at fixed positions in the sky, providing digital images every 10 s.If the same cloud detail could be identified in overlapping areas of stereo images taken with at least two cameras simultaneously, its three-dimensional position could be determined (Seiz et al., 2007).In our case, the horizontal distance between the cameras was too large for simultaneous observations, but cloud features moving with about constant speed across the camera view-fields could be used for photogrammetric analysis.The results are used to identify the causing aircraft, contrail positions vs. time with respect to NWP data, and to deduce contrail and atmospheric properties.For a direct comparison of simulated contrails with camera observations, we map the computed contrails on synthetic camera images.This requires algorithms which compute spherical coordinates for given pixel coordinates and their inverse.
The pixel positions of identifiable objects in digital camera images can be determined with standard image processing software.However, images of wide-angle cameras usually are distorted considerably (Weng et al., 1992;Garcia et al., 1997).In our case the distortion becomes obvious because straight contrails appear increasingly curved when coming close to the camera position, in particular near the edge of the camera image.The mathematical algorithm which describes the transformation of image coordinates into horizontal spherical coordinates or vice versa, including corrections for distortion is called a camera model.
Camera models have been used widely for astronomical observations, mainly for dark sky imaging, e.g., for meteor trace analysis (Oberst et al., 2004), also for gravity wave analysis in mesospheric airglow images (Garcia et al., 1997), noctilucent cloud observations (Baumgarten et al., 2009), cloud mapping using calibration with stars and aircraft with known positions (Seiz et al., 2007;Shields et al., 2013), or for automatic identification of stars in digital images (Klaus et al., 2004).
In contrast to stars, cloud features are more fuzzy and variable in time, and hence only allow less accurate geometric observations.Contrail features cover typically observation angles of one or a few degrees.Therefore, our camera model uses a simplified imaging geometry and distortion model and exploits symmetry assumptions to cover the full image frame.
This article describes camera models for two types of wide-angle cameras (a whole-sky imager and a more narrow one).The camera models correct for radial distortion inside the camera and for the orientation of the camera with respect to the horizontal coordinate system.The camera models are calibrated by using observations of the Sun, the Moon, planets and a few bright stars and landmarks.Moreover, we report results from aircraft track and contrail motion analysis with the camera models, and compare them with air traffic movement data, numerical weather prediction data, and simulations of contrail trajectories and spreading with a Lagrangian contrail model.

The cameras used
Cloud images were obtained in this study with four commercial digital video color cameras (Table 1).Photogrammetric analysis is described in detail for two of them (Figs. 1 and 2): 1.A wide-angle camera of type Mobotix D24M with L22 lens, installed on the roof of the Institute of Atmospheric Physics of the Deutsches Zentrum für Luftund Raumfahrt (DLR, German Aerospace Center) at Oberpfaffenhofen (OP).The wide-angle lens covers a limited field of view and faces roughly westward with some upward tilting.
2. A whole-sky fisheye camera of type Mobotix Q24M installed on the roof of the Meteorological Institute of the Ludwig-Maximilians-Universität in Munich (MIM), pointing vertically.
The local time of observation, received from an internet time server, is recorded with an accuracy of about 1 s.The two further cameras, MAY and HOP, are also of type 1, as the OP camera.

Camera model algorithms
The camera model provides the relationship between celestial azimuth and elevation angle coordinates (A, E) and pixel coordinates (X, Y ).X counts image columns from left to right, from 1 to 2048 and 1 to 1280, for camera 1 and 2, respectively; Y counts image rows from top to bottom, 1 to 1536 and 1 to 960.The azimuth A varies from 0 to 360 • (0 • : North; 90 • : East).The elevation is zero at the horizon and 90 • at the zenith.The algorithm makes use of the pixel and celestial coordinates of the camera midpoints, (X 0 , Y 0 , A 0 , E 0 ).For camera 1, the coordinates X 0 and Y 0 are set to the midpoint of the image plane and the angles A 0 and E 0 are found from astronomical observations near this midpoint.For camera 2, the midpoint is set to coincide with the vertical direction (E = 90 • , A = 0 • ) and the values of the coordinates X 0 and Y 0 are found by minimizing the model residuals, i.e., the root mean square (rms) differences between observed and computed star observations.A large focal length of the cameras is important for high resolution (Table 1), but its value is not used explicitly in the models.
The camera model describes the mapping between pixel coordinates (X, Y ) in the image and (X , Y ) coordinates in an imaginary plane onto which the spherical object coordinates (A, E) are projected (the so-called projection plane).The two camera types differ in their mapping transformations (see Fig. 3).For camera 1, which is tilted upward, the projection plane is the plane tangential to the celestial unit sphere, with the tangential point being at the image center.The rectangular coordinates in this plane are defined by gnomonic projection.For camera 2, the lens projects the sky onto a horizontal finite circular image resolved by a rectangular set of image pixels.Here, we choose the horizontal plane as projection plane, with the position angle of an object's projection point being its azimuth, and the center distance being proportional to its zenith distance angle.An affine transformation is constructed between the two Cartesian coordinate systems.Additionally, the camera model accounts for the distortion caused by the camera lens, i.e., its deviation from the ideal geometry of a so-called pinhole camera (van de Kamp, 1967).We assume that radial lens distortions are symmetric with respect to rotation around the image center.Instead of a polynomial function (e.g., Weng et al., 1992), we use an exponential function because of well-defined asymptotic behavior for small and large radius arguments.
With these definitions, we compute the parameters of the camera models as follows.First, we determine the parameters of the distortion function by comparing for all observations their center distances in the image with those in the projection plane.We then use the distortion function to rectify the observed pixel coordinates.Finally, following the classical Turner method (Turner, 1894), we construct an affine transformation between the rectified pixel coordinates and the computed locations of the celestial objects in the projection plane.This transformation represents the camera orientation and image scaling.Since the problem is over-determined, least-square fits are applied in the parameter computations for both the distortion function and the affine transformation.
Based on these considerations, the forward transformation is constructed as follows: 1.The camera pixels (X, Y ) are transformed into rectangular coordinates (x d , y d ) relative to the camera center, with y d pointing upwards in the image, (2) 2. The coordinates (x d , y d ) are mapped to (x , y ) assuming a radially symmetric image distortion, The ratio r/r d of the radii, with coefficients a > 0, b ≥ 0, c > 0, is used to correct for radial distortion, 4. For camera 1, the projection plane coordinates (X , Y ) are related to the angles (A, E) by trigonometry (van de Kamp, 1967): For camera 2, we use:  The inverse transformation is set up consistently as follows.
1. Camera 1 uses the gnomonic projection of the angles (A, E) to camera image coordinates (X , Y ) (van de Kamp, 1967): Camera 2 uses the inverse of Eq. ( 10), For both cameras, we compute the inverse of Eq. ( 5): The inverse solution r d (r) of Eq. ( 3) is determined by a Newton iteration, starting from a first guess r d = r/a.Finally the pixel coordinates are given by In total, both cameras use 11 free model parameters: Â, B, Ĉ, D, Ê, F , a, b, c, E 0 , A 0 for camera 1, and X 0 , Y 0 instead of E 0 , A 0 for camera 2. To avoid a conflict in scaling by a and the Turner coefficients, we normalize the latter such that the determinant Â Ê − B D = 1.This reduces the number of free parameters to 10.For their calibration, we assume that a set of observations (A i , E i , X i , Y i ), i = 1, 2, . . ., n, is available (see Fig. 3).In the following we will show how these parameters are determined.
For each observation, we compute R(A i , E i ) from Eqs. ( 2)-(6) and R(X i , Y i ) from Eqs. ( 12)-( 15).Ideally, both should be equal.Here, we determine the fitting parameters a, b, c entering Eq. ( 3) such that the differences have minimum rms values, Using Eqs. ( 3) and (4), we compute r i , x i , y i for each observation.Together with the X i , Y i from above, these data should satisfy Eq. ( 5): with unknowns Â, B, Ĉ, D, Ê, F .This set of equations can be expressed in matrix notation as For n > 3, Eq. ( 21) is over-determined.A least-square solution for the model parameters is found by solving the normal equations, These linear systems with 3 unknowns each are solved numerically.
Optimal values of A 0 , E 0 (camera 1) or X 0 , Y 0 (camera 2) are found when minimizing the properly weighted sums of the squared residuals of the observations for angles and pixel coordinates together with the radial residuals (Eq.19).The camera models are available as implementations in MS Excel 2010, Python, and Fortran, using various optimization algorithms (e.g., Schumann et al., 2012).

Camera calibration observations
For the determination of the camera model parameters (Table 2), we need a set of coordinates (X i , Y i ), i = 1, 2, . . ., n, related to known azimuth and elevation (A i , E i ).Basically, n = 3 linearly independent pairs of (precise) observations are sufficient to determine the 6 coefficients of the affine transformations (Eq.5).At least 3 (partly the same) observations spaced over the radial range are needed to determine the 3 coefficients in the radial distortion function (Eq.3).However, to compensate for observation errors, a large set of observations covering the whole field of view is desirable, n 3.
Here, we describe calibration for cameras OP and MIM.For these cameras, we use occasional sightings of celestial objects in images taken at times with clear sky.Since the orientation of the cameras is fixed, sightings taken at different, accurately measured times can be used for a single camera model.
In the actual application, the field coverage is far from optimal.One of the cameras was operated routinely only during daytime.Moreover, because of strong stray light at the urban observation positions, only a few usable nighttime observations are available.Therefore, the cameras are calibrated using observations of the Sun, Moon, and Vega for OP, and the Sun, Moon, Venus, Jupiter, and Sirius for MIM.
For these sightings, the coordinates of the brightest pixel are measured at high magnification using the software Irfanview.The planetarium software Guide 9.0 provides highprecision spherical coordinates with respect to the horizon, including the correction for atmospheric refraction, for the times when the digital images were taken.
For OP, the sightings are concentrated in essentially two image bands (see Fig. 3).One band (mainly Sun, some Moon and some Vega observations) runs from the upper-left to the lower-right corner and crosses the middle of the image.The other one (Moon observations) is parallel to the first and lies For MIM, observations are available mainly in the South.The northern range of azimuth angles between 292 • and 72 • is not covered.Three landmarks at low elevation angles at distances between 200 and 4500 m are used again for tests (see Fig. 4).
As a result, with coefficients as given in Table 2, the image resolution in degree per pixel, as derived from derivatives such as ∂A/∂X, is 0.052 in A and 0.045 in E at the image midpoint in OP, and 0.158 in E at the zenith and 0.109 and 0.215 in A and E at the horizon of the MIM camera.The resolution for MIM is slightly better than what was reported for multispectral whole-sky cameras before (Feister and Shields, 2005;Seiz et al., 2007).
The resolution is sufficient to observe contrails at = 100 m geometric scales up to a range radius R = N X 180 • /(π A) of more than 22 km (see Table 1 and Fig. 2).Here, N X is the number of image pixels in horizontal direction, and A is the azimuth angle range of 101 • , 66.3 • , 95.2 • for OP, MAY, HOP, respectively.For cameras OP, MAY, and HOP, this range is computed for elevation E 0 .For MIM, we compute the range at zero elevation with N X / A = ∂X/∂A = 9.14.

Discussion of the model accuracy
The model accuracy is limited by two factors: 1.The model is not a perfect description of the camera.
For example, the image distortion function approximates an a priori unknown relation between measured and real image center distances.A tangential distortion, e.g., (Weng et al., 1992), is not taken into account.
2. The astronomical sightings are affected by measuring errors in the images.The glare (or blooming effect, Seiz et al., 2007) caused by bright objects (Sun and Moon) and by the lens may cause errors of the order of several pixels, in particular close to the image borders.Uncertainties in the celestial coordinates provided by the planetarium software are several orders of magnitude smaller and can be neglected.
The model accuracy is assessed by testing how well the model matches the astronomical sightings it is based on.Since there are many more sightings than model parameters, residuals will be small only if the model is a good description of the real camera setup.The residuals of the camera model are computed from differences between given and computed image coordinates.Both the maximum and the rms residuals are evaluated (see Table 3).Further independent tests will be provided by sightings of aircraft with known positions (see Sect. 3.3).The rms angle residuals are in the range of 0.2 • or smaller, and the pixel residuals in the range of 3 (see Table 3 and Fig. 4).The pixel residuals for MIM are smaller and the angular residuals are larger than for OP because of coarser angular resolution per pixel.Given the short focal length of the cameras, and considering the fact that the angular diameter of the Sun and Moon is about 0.5 • , these residuals are within the expected measurement errors.Note, an angle error of 0.2 • corresponds to 35 m displacement at 10 km altitude in the zenith above the MIM camera.The pixel residuals are larger, mainly because of glare, but comparable to the accuracy of cloud feature observations.The landmarks are at very low elevation angles and at small distances, which could be a source of additional uncertainty, but the residuals are in the same range as for the other observations.Landmarks are time-independent and they show no systematic A residuals, which would arise from astronomical observations if the image time readings were systematically high or low.
To show the sensitivity to the distortion corrections (Eq.3), we applied the model also without the corrections.In this case, the rms residuals become about 20 times larger.The radial corrections mainly reduce elevation residuals, while the affine transformation mainly impacts azimuth values.The importance of the radial transformation can also be seen from the factor b (Table 2) which amounts to about 3.5 %.The exponential term becomes large for pixel radii larger than 1/c, of about 200 to 500.
While small residuals are a good indicator for the model accuracy within the image areas covered by observations, they cannot be used to assess the model quality in places where no observations are available.In those places the models are extrapolations based on symmetry assumptions.An important assumption is that the camera lenses exhibit a circular symmetry.For the radial distortion function an exponential function was chosen which is free of artificial oscillations.This would be different, e.g., if a polynomial had been used.Tests have shown that the difference of prediction and measurement (in pixels) for the image distortion function, Eq. (3), does not exceed the typical measuring errors of a few pixels independent of the position angles in the image.
Finally, the observational data have to be sufficient for a stable computation of the affine transformation (Eq.5).The parameters in this transformation become ill-defined if the image area covered by observations degrades to a line.Fortunately, for both cameras the covered areas have large extensions in both coordinate directions.For modern cameras, we may expect nearly non-skewed and isotropic pixel orientations, so that D ≈ − B, Ê ≈ Â.Here, Â and Ê differ by less than 0.01 % for both cameras.Small values of B and D are expected if the camera mounting is precisely aligned horizontally.These relationships can also be used to assess the quality of the model input for a limited set of observations.

Transformations between observation angles and geographic coordinates
For relating geographic positions of an object to camera observation angles, we use Cartesian coordinates (x, y, z), in m, with the horizontal plane x-y tangential to a sphere, approximating the Earth with mean radius R ≈ 6371 km, at longitude λ = λ C , latitude φ = φ C , and altitude z = z C of the camera C above mean sea level (a.s.l.).Here, x, y are the orthogonal horizontal geographic coordinates in eastern and northern directions and z is the vertical coordinate a.s.l.For small distances λ − λ C and φ − φ C , the coordinates x, y are related to λ and φ, in degrees, approximately by For large λ − λ C and φ − φ C , great circle computations (van de Kamp, 1967;Earle, 2005) are required instead.
In the computation of the elevation angle E and the projected distance d on the ground between the object and the camera we account for the curvature of the Earth surface.Tests have shown that the curvature may be ignored for contrail altitude and wind speed determinations for altitudes above 8 km and distances below 50 km.For larger distances and low elevation angles, however, the curvature must be considered.An object at altitude H sinks below the horizon (E = 0) at the horizontal distance a = √ 2RH + H 2 ; e.g., a = 35.7,112.9, 357.1 km for H = 0.1, 1, 10 km, respectively.
For an object P at altitude H +z c a.s.l.viewed from camera C at angles A, E above Earth horizon with Earth origin in O and effective Earth radius R = R + z C (see Fig. 5), the general triangle OCP is defined by two given side lengths OC and OP and one angle, ψ = d/R or α = E + 90 • .Hence, we find the distance d along the Earth surface at sea level and the geographic coordinates (x, y) = g(A, E, H ) from Similar relationships to compute A, E for given x, y, H are given in Garcia et al. (1997), but ours are simpler and allow for explicit inversion to compute x, y for given A, E, H .

Observations
We observed four crossing contrails passing the view field of camera OP at about 08:42 UTC 3 November 2012 (see Fig. 6) (all the clock times refer to UTC).The air was clear with at least 100 km visibility.A westward wind was strong, with about 50 m s −1 (e.g., according to ECMWF data) at the contrails' altitudes, so that the contrails happened to move into the direction of Munich.Because of the cross pattern, the same set of contrails was clearly identified at MIM about 10 min later.The contrails are named C1, C2, C3 and C4, according to their clockwise appearance in the OP images.These observations will be used to determine the contrail altitudes, tracking speeds, and widths.In addition, a short-lived contrail C5 was spotted.These observations also provide information on the humidity.
The principle of this analysis can be understood from Fig. 7.The contrail (here C1) is observed by the two cameras at various times at various elevation angles between about 30 and 90 • , approaching OP from the west and passing MIM over a distance of about 50 km.The mean altitude and the mean speed of the contrail along the x axis are determined by a fit to the observed elevation angle changes with time.Details are described in Sect.3.2.
In addition to the contrails, some of the related aircraft were visible in the camera images.For contrails C1, C3, and C4, the contrail-causing aircraft could be detected in earlier images, either by spotting the aircraft themselves or their fresh trailing contrails.The times of first visibility and related information are listed in Table 4.The aircraft causing contrail C2 was not visible, but the first detection of C2 west of OP at about 08:30 could be traced backwards with wind to identify the aircraft that caused this contrail in air traffic data a few minutes earlier.
The contrails were visible in the MIM images until about 09:09, i.e., for about 40 min.During this time, the contrails grew in width and got advected with the winds over a distance of about 120 km.The contrails appeared to become optically thicker with time (measurements of the optical depth would require multispectral cameras, Seiz et al., 2007).Hence, all these contrails are classified to be persistent.
The same contrails were incidentally observed by a high resolution camera MAY from Munich (no. 3 in Table 1), 1.27 km north of MIM.This type 1 camera has a rather narrow field of view with less distortion.Camera model 1 provides a reasonable approximation for this camera even without corrections for radial or linear distortion.Since this camera is not in fixed position, we estimated orientation and scaling from 3 landmark and 3 Sun observations.The camera pointing accuracy is estimated to about 1 • , as supported by the observation of an aircraft with a shortly visible contrail, the position of which is given by ADSB data.
Finally a low-resolution webcam HOP (no. 4 in Table 1) of the Observatory Hohenpeissenberg (Deutscher Wetterdienst), 37.53 km southwest of OP, calibrated with a few Sun, Moon, star (Arcturus and Vega), and landmark observations, documents the scene in westward direction.

Contrail altitude and wind speed
Without knowing the aircraft waypoints, the contrail altitude z can be determined by a least-square fit assuming the observed contrail positions can be tracked along lines at constant altitude with constant wind speeds U, V in east/north (x, y) directions.For this purpose, we manually measure coordinate pairs (X, Y ) for selected points within the contrail in images taken at various times and locations (here in 3 images from OP and 3 from MIM at times between 08:39 and 08:56 UTC).For wider contrails we use the higher (in elevation for OP, or most eastward for MIM) contrail edge to capture the contrail tops.From the pixel coordinates, the azimuth and elevation (A, E) pairs are computed using the camera model (Eq.1).With estimated altitude z (e.g., 10 km), the distance d and the horizontal geographic coordinates (x, y) are computed from Eq. ( 26), and then the altitude is corrected iteratively to minimize the rms to the observations, as described below.The final results are shown in Fig. 8.The contrail coordinates (x, y) are assumed to be close to straight lines The values C and s at various times are fitted so that the sum of (x−x c ) 2 +(y−y c ) 2 over all points along a contrail in a single image assumes its minimum.The plots in Fig. 8 show that the contrails are in fact very close to linear.This also shows that the camera model correctly removes the distortion of straight lines.For a moving line without marked features, one can determine the tracking speed only in the normal direction.Here, we follow the advection of the cross point between the contrail lines and and the line connecting the two camera positions with slope S = (y MIM −y OP )/(x MIM −x OP ).For given nonzero slope difference s − S, this point moves with constant advection speed U a = (s U − V )/(s − S).Note that the computed geographic coordinates X, Y depend on altitude H = z − z C by Eq. ( 26).The altitude z a.s.l., the speed U a , and the reference coordinate x 0 are free parameters which are determined with an optimization routine such that x c (t) = x 0 + U a t agrees with the observed contrail line crossing point in the least square sense.This works for slopes s significantly different from S. In this case, we have three unknowns (z, x 0 , U a ) and six measurements (x c at the six observation times), hence this minimization process has a well-defined solution.For contrails nearly parallel to the he lines show the linear fits to the contrail positions as computed from the observed image pixles (X, Y ) using e camera models (A, E) = F (X, Y ) and the geographic coordinates (x, y) = g(A, E, H) for fitted geometric ltitudes z = zC + H a.s.l.Grey lines show lake positions for orientation.

Fig. 7.
Viewing directions to contrail C1 from cameras OP and MIM at a sequence of times (08:39:08, 42:08, 42:38, 48:42, 51:43, 55:15 UTC, from left to right).The contrail altitude is the result of the fit described in Sect.3.2.The contrail x values are the positions where the contrail line cuts the x axis through OP. connection line between the cameras (C4 in our example), we need to identify contrail features (here the end points) which can be assumed to move with constant wind speed (U, V ).
Table 4 lists the results.The fit results are accurate of up to 230 m rms errors for z, and 0.23 m s −1 for advection or wind speed.This was found out by systematically repeating the analysis with random selections of subsets of the camera readings X, Y .The rms error is largest for C3, which has an orientation not far different from the wind direction.
The contrail altitudes derived from the camera fits can be compared with the aircraft flight level information (see Table 4).Note that the camera observes the geometric altitude.Aircraft flight levels (FL) are pressure altitudes z p in hft defined for static pressure in the International Civil Aviation Organization (ICAO) standard atmosphere (ICAO, 1964).In our case, with lower-than-average surface pressure and warmer atmosphere up to 8 km, the geometric (or geopotential) altitude z g is lower than the ICAO pressure altitude z p , according to ECMWF data, for this case ( z ICAO =   OP and MIM (times 08:39:08,42:08,42:38,48:42,51:43,55:15 UTC) derived from the observations in images from cameras at OP and MIM.The lines show the linear fits to the contrail positions as computed from the observed image pixels (X, Y ) using the camera models (A, E) = F (X, Y ) and the geographic coordinates (x, y) = g(A, E, H ) for fitted geometric altitudes z = z C + H a.s.l.Grey lines show lake positions for orientation.z g − z p = −200 ± 30 m).Hence, the effective altitude difference is z g − z p − z ICAO .With this correction, the altitude differences are within ±100 m (see Table 4).The wind speed components U and V derived for case C4 at 8.7 km altitude agree within 5 % with the ECMWF data (44.5 and 6 m s −1 ), see Fig. 13.
Errors of the order of 200 m my be acceptable when considering other sources of uncertainty: The lower part of a contrail sinks during the wake vortex phase by a range of the order of 50-300 m depending on aircraft and atmosphere parameters (e.g., Schumann et al., 2013).Differences may also result from uncertain pixel readings for thick contrails, horizontal variations in the wind speed (the U wind component seems to increase with y), and atmospheric wave motions.However, the accuracy of the derived contrail altitudes is consistent with stereo camera cloud altitude results (Seiz et al., 2007).
For analysis of contrail widths, we use overlays of horizontal x-y grids into the image, as shown for example in Fig. 9.The geographic grid coordinates x, y, z are specified for given contrail altitudes and given horizontal resolution.The view angles A, E and the image coordinates X, Y are computed using Eqs.( 1) and ( 29).The widths observed over OP and MIM, as listed in Table 4, are determined by matching the observed contrails with the grid, with about 100 m accuracy.The width of the short-lived contrail C5 is at the limit of resolution.For the others, the width accuracy is limited mainly by the contrail shape and contrail edge contrast against clear sky.With time, the persistent contrails become wider.The width for C3 includes the sum of primary and secondary wake parts which can be visually distinguished in this case.Because of positive wind shear, the primary wake appears at the more westward edge.

Aircraft identification
For contrails C1, C3, and C5, the pixel coordinates X, Y of the aircraft sightings were measured (typically with ±2 pixel uncertainties).For C4, the first contrail appearance was located in the OP camera images.These data were converted into azimuth and elevation angles A, E using the OP camera model, Eq. ( 1).With this information and an estimated altitude z, we use Eq. ( 26) to estimate the geographic horizontal coordinates (x, y, z, t) of the first contrail sightings.
From the German air traffic control agency (DFS, Deutsche Flugsicherung), we obtained the waypoint coordinates of all aircraft movements above about 7 km for this day over Germany.The data give the waypoint coordinates (x, y, z, t) in 1 min intervals.The DFS positions were compared with ADSB observations (available every 5 s).Presently, most aircraft in operation are equipped with ADSB transponders.Exceptions may occur in particular for small jets.The ADSB data cover all the flights for which contrails were identified in this case study, and give (within round-off or time interpolation errors) identical position values.For the following analysis, DFS data are used.
Plotting the coordinates x, y, z, t of the first contrail sighting together with the coordinates of the aircraft waypoints, the aircraft flights could be identified without doubt even for  rough altitude estimates (about 10 km).With altitude from these data, the geographical coordinates of the contrail sightings were matched accurately.For example, Fig. 10 shows the positions of aircraft sightings together with the flight coordinates in an x-y plane.For given altitude, the horizontal aircraft positions as derived from the camera observations and from the waypoint data agree within 200 m, or better than 1 %, even at the most remote distances.This demonstrates nicely the accuracy of the camera model.
For C2, the aircraft was too far away (more than 70 km) to be visible in the photos.Here, Fig. 10 depicts the position of the first contrail sighting.An aircraft flying further west about 4 min before the first sighting caused contrail C2.
The position of the first appearance of C4 agrees accurately with the aircraft track (in horizontal position, time and altitude).The aircraft causing C4 was climbing while flying westward.The aircraft flight level listed below is that at the time of first contrail appearance.
An aircraft with the short contrail C5 (about 1 to 2 km length, i.e., less than 10 s maximum age) is visible in Fig. 6, at least in the full-resolution original images.From the observed coordinates and the waypoint data, the aircraft was clearly identified (see Fig. 10).

Synthetic contrail images
For given aircraft waypoint and wind information, the trajectories of the contrail waypoints can be computed using the Lagrangian advection part of CoCiP (Schumann, 2012) (see Fig. 11).
For each aircraft waypoint (x 0 , y 0 , z 0 , t 0 ) along the flight track, the local wind components are linearly interpolated from the NWP data in time and space.The NWP data   age mainly because of differences between NWP-derived and true wind speeds.Note, that 1 m s −1 wind error for a contrail age of 1000 s implies a position error of 1000 m, which corresponds to about 6 • angular displacement for an overhead contrail at 10 km altitude.If contrail C1 were computed for a 15 s later time, its position would agree perfectly with the observation in the MIM photo.It seems that the true wind was slightly stronger than predicted by the NWP model, both in x and y directions.
Figure 11 also depicts the left and right boundaries of the contrails by plotting two lines at the same altitude as the contrail center line, shifted laterally by the half widths in geographic space.The computed and observed widths agree fairly well for C1, C4 and C5, but the observed contrails C2 and C3 are about a factor of two wider than simulated, possibly because of underestimate of small-scale shear in the NWP data.Contrail C4 experiences the weakest shear and may stays more narrow, therefore.Hence, such synthetic images open a new approach to test and possibly improve contrail modeling.
Synthetic contrails C1 to C4 were plotted also for the two other cameras.From HOP one of the four contrails (C2) was visible (at about 08:27) in reasonable agreement with synthetic images.From MAY, the contrail cross was observed and simulated while passing towards East (see Fig. 12).The picture supports the approximate validity of the synthetic contrail positions, and the persistence and increasing width of the contrails, besides many other interesting cirrus structures.

Checks of humidity data and Schmidt-Appleman threshold
Contrail and cirrus properties are strongly sensitive to relative humidity.Observations and numerical humidity predictions are difficult for many reasons.The formation of icesupersaturation depends on vertical motion, temperature and cirrus ice microphysics (Tompkins et al., 2007).Layers of ice-supersaturation are often rather thin (Gierens et al., 2012) and hence difficult to resolve numerically.
Figure 13 shows wind, relative humidity over ice (RHi), and temperature vs. altitude as computed from ECMWF data.The model predicts ice-supersaturation between about 9.0 and 11.3 km pressure altitude, with a local minimum in RHi near 10.3 km altitude.For the NWP values of temperature, humidity and pressure, and for aircraft burning kerosene with an overall propulsion efficiency of 0.3, the Schmidt-Appleman criterion (SAC) implies contrail formation for pressure altitudes z in the altitude range 9.5-16.5 km.The contrails C1 to C5, with the exception of C4, formed in this altitude range.C4 formed at about 0.7 km lower altitude.This indicates a possibly higher ambient humidity at this level than predicted by ECMWF.
At the altitudes of the observed persistent contrails, C1 to C4, the RHi must have been above 1.This is indicated by the    (Schumann, 1996).η measures the work performed by the aircraft engines by thrust and true air speed for given combustion heat and fuel flow per time unit.For cruising jet aircraft, η is typically between 0.3 and 0.38. Figure 13 shows T SAC,1 , computed from ECMWF values for pressure and RHi and for η = 0.3.
In this case, contrail C4 could not be explained.The ambient temperature was about −39 • C, more than 7 K above the SAC temperature (−46.4 • C).The temperature accuracy of such NWP models is typically within 1 K (confirmed by comparison to other NWP output for this case).An increase of η by 0.1 corresponds to an increase in RHi by 33 %; both cause 1.55 K higher T SAC .Hence, even η = 0.4 would not suffice to make T SAC larger than T .During climb, as in this case, η is usually smaller than at cruise because of lower aircraft speed.Hence, the ambient humidity must have been strongly ice-supersaturated.In clear air, humidity may reach or slightly exceed the homogeneous freezing limit (Koop et al., 2000), which equals liquid saturation near T = −40 • C (about 1.45).Only with such high humidity, as indicated for C4 in Fig. 13, the atmosphere was just cold enough to let contrail C4 form as an exhaust contrail according to the Schmidt-Appleman criterion.
The short contrail C5 at 12.1 km indicates subsaturation at this pressure altitude.Hence, the NWP-predicted subsaturation at z > 12 km (see Fig. 13) is confirmed by this contrail observation.From the results for C1-C4, the layer with icesupersaturation reached over 8.6-11.8km pressure altitude, nearly 40 % larger than predicted (9.0-11.3km).

Conclusions
This paper describes methods for contrail tracking and analysis of contrail properties from video camera observations, in particular contrail geometric altitudes, widths, and motion speeds.The methods are applied to a case study of contrail observations using two different kinds of wide-angle video cameras (whole-sky imager with fisheye lens or wide-angle cameras with smaller field of view), placed several kilometers apart.
Photogrammetric methods are described for the two camera types.The camera models allow us to determine azimuth and elevation angles for given image coordinates and vice versa.The models account for linear and radial distortions.For the calibration we use mainly sightings of bright celestial objects together with some landmarks and aircraft with fresh contrail observations.An incomplete coverage of the field of view with such observations is overcome by exploiting reasonable symmetry assumptions in the camera models.The accuracy of the models, demonstrated by the residuals between analyzed and observed coordinates and by the agreement of the observed positions of young contrails with waypoint data, is within the range of the expected measuring errors.
The case study describes a "4-contrail cross" persisting for about 40 min, together with a short-lived one.Some of the contrail forming aircraft were visible and identified by comparison to air traffic waypoint data.The waypoint information from DFS and from ADSB data was found to be in good agreement.The other contrails were related to aircraft flight tracks by means of contrail trajectories.From the comparison of observed positions with movement data, we found that the camera models and observations with two cameras allow determining the altitude and horizontal position of the contrails to an accuracy of better than 230 m, width to about 100 m, and the mean horizontal tracking speed to about 0.2 m s −1 .In comparing altitudes, differences between ICAO standard atmosphere pressure altitudes and geometric altitudes are significant.
The observed contrail evolution is compared with simulated contrails.Contrails are simulated with the contrail prediction model CoCiP, a Lagrangian model using air traffic movement data and numerical weather prediction (NWP) data as input.The results are projected on camera images.Here, the availability of the inverse camera model was essential.
The presence of a contrail constrains the relative humidity being below or above the thresholds required for contrail formation (Schmidt-Appleman criterion) and contrail persistence (ice-supersaturation).The observations show spreading contrails, apparently with increasing optical depth, suggesting ice-supersaturated ambient air at contrail pressure altitudes (from 8.7 to 11.7 km).The ice-supersaturated layer is found considerably thicker than predicted by the NWP model used.In fact, to understand contrail C4 as being formed as an exhaust contrail, the aircraft must have flown in air with high relative humidity, close to liquid saturation at the time of contrail formation.
The model tends to underestimate the contrail widths, indicating underestimates in the initial contrail depth or ambient shear (from the NWP data) and turbulent mixing.With age, the horizontal contrail positions become increasingly sensitive to the assumed wind field.Although the camera derived wind data agree with ECMWF data within about 2 m s −1 , such small differences cause notable shifts of aged contrails in the camera images.
Hence, multiple camera observations of contrails provide insight into contrail dynamics and tests of numerical weather prediction and contrail models.Further insight into cloud dynamics and NWP may be obtained by using these and similar cameras for longer time periods and possibly at additional sites.In the future, one may envisage using a dense network of fixed-mounted video cameras, preferably in connection with other remote sensing methods, to observe contrails and wind fields, and to determine constraints for the humidity over a larger region and longer time period.Even single camera observations may be useful in this respect when combined with altitude information from other sources (e.g., ADSB data or ground-based lidar).The work described in this paper was initiated by observations of a special cirrus cloud which looked similar to contrails, but was not easily attributable to specific aircraft flights.The analysis of this special cirrus with the given photogrammetric methods and further observations is to be described in a future paper.

Fig. 1 .
Fig. 1.Wide-angle camera in Oberpfaffenhofen (OP, left panel) and Munich (MIM, right panel).The sphere in the top-cross of the St. Markus church, above the MIM camera, is used as southeastern landmark.

Fig. 2 .Fig. 1 .
Fig. 2. Positions, view angles and view ranges of 100 m resolution for the cameras Oberpfaffenhofen (OP), Munich (MIM), private (MAY), and Hohenpeissenberg (HOP) in Southern Germany.MIM and MAY are 1.27 km apart, so nearly coincide in this figure.The region is located between Lake Constance (West), Chiemsee (East), Danube (North), and Inn (South).The lakes Ammersee and Starnberger See are within the range of MIM, south of Munich.

x
= x d r/r d , y = y d r/r d .(4) 3. The rectified image coordinates (x , y ) are mapped with an affine transformation to projection plane coordinates (X , Y ) which accounts (as discussed in Sect.2.3) for camera inclinations and rotations (parameters B, and D), horizontal shifts (parameters Ĉ, and F ) and scaling (parameters Â, and Ê), possibly different in the two directions, X = Âx + By + Ĉ, Y = Dx + Êy + F, )5.Finally, equations sin(A ) = S A , cos(A ) = C A , with given S A , C A , imply if C A > 0: A = sin −1 (S A ), else: A = 180 • − sin −1 (S A ).Negative values of A are incremented by 360 • .the top-cross of the St. Markus church, above the MIM camera, is used as southeastern landmark.

Fig. 2 .Fig. 2 .
Fig. 2. Positions, view angles and view ranges of 100 m resolution for the cameras Oberpfaffenhofen (OP), nich (MIM), private (MAY), and Hohenpeissenberg (HOP) in Southern Germany.MIM and MAY are 1.27 apart, so nearly coincide in this figure.The region is located between Lake Constance (West), Chiemsee (E Danube (North), and Inn (South).The lakes Ammersee and Starnberger See are within the range of MIM, s of Munich.

Fig. 3 .Fig. 4 .Fig. 3 .
Fig. 3. Coordinate lines of mutual transformation for the two camera types (left panels: limited wide-angle, right panels: whole-sky imager), and field coverage with astronomical sightings (circles) for camera type 1 (e.g.OP, left panels) and type 2 (MIM, right panesl) in image coordinates (X, Y , top panels) and celestial angle coordinates (A, E, bottom panels).The grey triangles denote landmark observations.

U.
Schumann et al.: Contrail camera observations closer to the lower-left corner.Measurements in the upperright or lower-left image corner are missing.Observations of landmarks at low elevation angles at distances between 130 and 360 m are used for independent model tests (see Fig. 4).

Fig. 5 .
Fig. 5. Illustration of angles E, β, ψ, distance d, altitude H and effective Earth radius R = R + z C for an object at position P and camera at C.

Fig. 5 .
Fig. 5. Illustration of angles E, β, ψ, distance d, altitude H and effective Earth radius R ′ = R + ∆z C for an object at position P and camera at C.

Fig. 6 .
Fig. 6.A "4-contrail cross" formed by contrails C1, C2, C3, and C4 west of Oberpfaffenhofen and over Munich between 08:39:08 and 08:55:15 UTC 3 November 2012.C5 is a short-lived contrail.Contrails and other clouds appear white in front of the blue sky.The camera in OP is oriented westward covering elevations of about 0 to 60 • .The fisheye camera in MIM is oriented towards the zenith, with North, East, South, and West at the top, right, bottom and left image boundaries, respectively.To compare cloud patterns in the photos from OP and MIM, one has to rotate the MIM photos counterclockwise by about 90 • and to mirror along the West-East axis.

Fig. 9 .
Fig. 9. Section of Fig. 6 (red color part), from OP at 08:40:28 UTC, with the 5 contrails C1 to C5, w grids of geographic x-y lines for orientation, a coarse one (1 km grid spacing) at z = 11 km altitude a.s a fine one, rotated into flight direction of C5, with 100 m grid spacing, at z = 12 km.The insert in th left corner shows contrail C5 and the fine grid enlarged.

Fig. 10 .Fig. 8 .
Fig. 10.Aircraft flight tracks (red lines) in horizontal coordinates (x, y) relative to the positions of the in Oberpfaffenhofen (OP, cross) as identified from aircraft waypoint data.Sightings of the aircraft (black connected with black lines), based on pixel coordinates in the camera images and traffic data flight al The single black circles for C2 and C4 denote the first contrail sighting.The blue lines (linear fits of ob positions) locate the contrails C1 to C4 at times 08:42:38 and 08:51:43 UTC.Grey lines locate the lakes

Fig. 9 .
Fig. 9. Section of Fig. 6 (red color part), from OP at 08:40:28 UTC, with the 5 contrails C1 to C5, with grids of geographic x-y lines for orientation, a coarse one (1 km grid spacing) at z = 11 km altitude a.s.l., a fine one, rotated into flight direction of C5, with 100 m grid spacing, at z = 12 km.The insert in the lo left corner shows contrail C5 and the fine grid enlarged.

Fig. 10 .Fig. 9 .
Fig. 10.Aircraft flight tracks (red lines) in horizontal coordinates (x, y) relative to the positions of the cam in Oberpfaffenhofen (OP, cross) as identified from aircraft waypoint data.Sightings of the aircraft (black cir connected with black lines), based on pixel coordinates in the camera images and traffic data flight altitu The single black circles for C2 and C4 denote the first contrail sighting.The blue lines (linear fits of obser positions) locate the contrails C1 to C4 at times 08:42:38 and 08:51:43 UTC.Grey lines locate the lakes.

.
Contrail horizontal coordinates (x, y; symbols) for contrails C1 to C4 atOP and MIM (times 08:39:08,   , 42:38, 48:42, 51:43, 55:15 UTC) derived from the observations in images from cameras at OP and MIM.nes show the linear fits to the contrail positions as computed from the observed image pixles (X, Y ) using mera models (A, E) = F (X, Y ) and the geographic coordinates (x, y) = E, H) for fitted geometric es z = zC + H a.s.l.Grey lines show lake positions for orientation. 27

Fig. 10 .
Fig. 10.Aircraft flight tracks (red lines) in horizontal coordinates (x, y) relative to the positions of the camera in Oberpfaffenhofen (OP, cross) as identified from aircraft waypoint data.Sightings of the aircraft (black circles connected with black lines), based on pixel coordinates in the camera images and traffic data flight altitudes.The single black circles for C2 and C4 denote the first contrail sighting.The blue lines (linear fits of observed positions) locate the contrails C1 to C4 at times 08:42:38 and 08:51:43 UTC.Grey lines locate the lakes.

Table 1 .
The cameras.

Table 3 .
Maximum and rms residuals in angles A, E (in degree) and image coordinates X, Y (in pixels) and number of observations n.

Table 4 .
Observed and computed contrail properties with FL pressure-altitude and observed geometric altitude a.s.l.

Schumann et al.: Contrail camera observations a
SAC,1 is the Schmidt-Appleman criterion threshold temperature for contrail formation for the RHi obtained from ECMWF and an overall propulsion efficiency η = 0.3; T SAC,2 is the same for higher η (0.4) and RHi equal to RHi max .function of ambient pressure, fuel properties (combustion heat and water emission index, Q = 43.2MJ kg −1 , EI H 2 O = 1.23), and overall propulsion efficiency η • E, 48 • N), vs. pressure altitude a.s.l. at 09:00 UTC.The symbols in (a) denote the wind speeds derived from the camera observations for C4, the circles in (b) symbolize estimated RHi values for observed persistent or short contrails (grey or open) for the five contrails.In (c), T