Unbiased estimations of atmosphere vortices: the Saturn's storm by Cassini VIMS-V as case study

The size determination of dynamical structures from spectral images poses the question where to fix the shape’s boundary. Here we propose a method, suitable for nearly elliptical shape, based on the fit of a 2D Gaussian to the pixel intensities of the spectral image. This method has been tested on a vortex structure embedded in the wake of the 2010 Saturn’s giant storm. On January 4th 2012 the Visual and Infrared Mapping Spectrometer (VIMS), onboard Cassini, observed a giant vortex in the Saturn’s northern hemisphere. The structure was embedded in the wake storm system detected on December 2010 by Fletcher et al. (2011). Therefore, all the VIMS observations focused on the Saturn’s storm have been analyzed to investigate its morphology and development. VIMS detected the vortex from May 2011 up to January 2012. The shape and size evolution has been determined for the vortex cloud top, visible at 890 nm. The largest size resulted 4000 km about and seemed to shrinks continuously up to January 2012, while the shape varied in the second half of the year. The vortex oscillated in 2 degrees latitude around 37oN planetocentric latitude, and drifted in longitude by ∼ 0.75 deg/day in westward direction.


Introduction
Sunlight passing through the atmosphere is scattered and absorbed by gases and aerosols to produce the characteristically nebulous appearance. The amount of scattering and absorption that light undergoes depends upon the abundance of aerosols as well as the wavelength of the radiation. Light at certain frequencies interacts with the gaseous molecules and aerosols more effectively than at others; consequently, the depth to which light effectively penetrates and escapes the atmosphere is a function of wavelength. By using this wavelength dependence it is possible to image clouds at different altitudes. Clouds, depending from both their cold/hot nature and their growing mechanism, can assume a variety of shapes and sizes. In this work we propose a method suitable to estimate clouds quite elliptically shaped. Though this technique could be applied to any planetary cloud system, it has been setup aiming to infer the evolution of a particular vortex observed on Saturn in occasion of a giant storm sprouted at the end of 2010.
Ground-based telescopes detected the first signs of the storm on 5 December 2010 at 21:05:22 Universal Time (UT), but it was possible also to observe and probe it by the Cassini mission. Cassini Composite Infrared Spectrometer (CIRS) (Flasar et al., 2004), Radio Plasma Wave Science (RPWS) (Gurnett et al., 2002), Imaging Science Subsystem (ISS) (Porco et al., 2004) and Visual and Infrared Mapping Spectrometer (VIMS) (Brown et al., 2004) sounded different pressure levels and exploited a variety of features associated to the storm development. Many works have already been carried out on this topic by means of the CIRS thermal infrared observations (Fletcher et al., 2011(Fletcher et al., , 2012, the RPWS radio pulses (Fischer et al., 2011) and the camera images analysis (Sayanagi et al., 2013) in the 2011 January to October lapse time.
The storm structure, in upper atmosphere, can be usefully analyzed by ISS and VIMS. The instruments cover the UV (only ISS), the visible (VIS) and the Near InfraRed (NIR) spectral ranges. They observed the giant storm on Saturn during its whole life cycle, with some regularity. The first detections of the outbursts' tail started on February 2011, but the feature object of this report has been observed, in various occasions and from different distances, until December 2012. Here we will describe the position and morphology history of that feature from its first occurrence up to January 2012. We assume that the vortex monitored by VIMS is the same described in Sayanagi et al. (2013), for its position, similar shape and remarkable dimensions, then we analyze the series of measurements reported in Table 1, starting with the ISS observation on 12 January 2011 and ending with the VIMS observation on 4 January 2012. In section 2 we will describe the data selected for the analysis and the processing algorithms that produced the required outcomes. In section 3 we will discuss our results and summarize our work.

Observations and data handling
In this paper the observations by the VIMS visual channel and by ISS-WAC (Wide Angle Camera), listed in Table 1, are analyzed. VIMS is a hyperspectral imaging spectrometer; its visual channel (hereafter VIMS-V) has 96 bands ranging from 350 to 1046 nm with a 7 nm nominal bandwidth. ISS-WAC is equipped with two filter wheels with a wide variety of filters. Both VIMS-V and ISS-WAC raw data are publicly available from the NASA Planetary Data System (PDS) Imaging Node. Our search on VIMS-V datasets has been carried out principally on the archive of IAPS-INAF, where updated geometrically and radiometrically calibrated data are available for the Italian scientific team. ISS observations were selected with the purposes of filling the large gaps in time in which VIMS-V was used to acquire data for other scientific objectives, or when Cassini's orbit was not adequate to observe the relevant latitude. In Table 1 the list of the VIMS-V and ISS-WAC observations is reported together with the respective pixel sizes, obtained through the conversion of the images from a geographical to a Universal Transverse Mercator (UTM) projection. The projection calculations have been done exploiting some tools and facilities from ENVI-Exelisvis (www.exelisvis.com) and ESRI (www.esri.com) software. UTM projection provides a north-south oriented swath of little distortion. Each of these swaths is called a UTM zone and is 6 • wide in longitude. The first zone conventionally begins at the International Date Line (180 • , using the geographic coordinate system). ESRI's ellipsoid and datum approximations of the Saturn's irregular shape and the mathematical model have been used to describe the location of unknown points on Saturn's surface, respectively. This geographic information, along with the pixel geo-referencing by C-SPICE kernels (Acton, 1996), permitted us to evaluate the spatial frames' dimension and orientation. The raw VIMS-V data are radiometrically calibrated by using the official calibration pipeline (McCord et al., 2004;Filacchione et al., 2007). The VIMS and ISS data have been geo-referenced using ad hoc algorithms based on NAIF-SPICE tools. No radiometric calibration has been applied to ISS data since in this paper we only focused on the morphology (shape and size) evolution of the vortex and on its displacement in the examined time period. The ISS images inspection aimed principally to maximize the probability of a univocal recognition for the structure under analysis, adopting the criterion used by del Ro-Gaztelurrutia et al. (2010). We included two ISS-WAC measurements in our survey: on 12 January, a time period without VIMS-V observations, and on 24 August, when both instruments, in about half an hour and with the same observation geometry approximately, observed the same region. The 24 August observations have been used to compare the intensity scales of the two instruments and to check their correspondence. The two images have been resized on a geographic base, re-sampled to a common dimension of pixel number and normalized on a monochromatic color range scale (8 bit, 0-255 levels). In the resulting images ( Figure 1)

Results and conclusions
The evaluation of the limits of the feature under study presents the problem of identifying its boundary without ambiguities. Not all the spectral images at different wavelengths show the same appearance of the vortex. As can be seen in Figure 2, where three spectral frames of the January 2012 VIMS-V observation are shown, there is a great difference between the vortex shape as revealed at the ultraviolet wavelength (left panel, 373 nm), in correspondence of a methane window (central panel, 752 nm) and in the deepest methane band of the visible range (right panel, 890 nm).
Each one of these wavelengths sounds the atmosphere at a different depth. The UV wavelength (left panel, 400 nm) senses relatively deeply, but the color scale (corresponding to spectral radiance on the sensor) is influenced mostly from particles absorption and Rayleigh scattering by gases.
The 752 nm (central panel, high continuum reflectivity) sounds the deepest level and is the least affected by gas scattering. In the methane absorption band (right panel, 890 nm) the reflected light comes principally from the highest particles, because methane absorbs both the transmitted and the diffused radiation deepen in atmosphere. Thus the bright oval feature in Figure 2c corresponds to a high altitude light reflection, well-contrasted against a darker, deeper background caused by methane absorption.
This study is carried out on the third case because it offers a simpler structure showing the top vortex as a bright oval (BO in the following). However a quantitative criterion to choose where inserting the boundary line between the BO shape and the background had to be fixed also in this case. The VIMS-V and the ISS-WAC data have been normalized in order to make both datasets homogeneous by: where S is the signal (spectral radiance for VIMS-V and Digital Number (DN) for ISS-WAC), S m and S M the minimum and maximum values for each dataset, and S N is the normalized value, respectively. In Figure 3 the spectral frames at 890 nm for the observations listed in Table 1 before (top) and after (bottom) the re-sampling and normalization are reported. The nearly elliptical shape of BO at 890 nm has then been fitted with a two-dimensional, elliptical Gauss function in order to objectively determine its boundaries, sizes and shapes. The fitting equation F(x, y) = A0 + A1*exp (-U/2) has been applied to the observations reported in Table 1, but for the ISS W1692857065, referring to 752 nm. In the fitting equation U is the elliptical function and A0 and A1 are the intercept and the scale factor of the elliptical Gauss function. In the U equation a and b are the FWHM of the Gaussian in the X and Y direction (sample and line in the pixel coordinate system), respectively. A code from the IDL library has been used for this calculation. This procedure requires operating on a rectilinearly gridded data ( Figure 4). The fitting routine is based on an adaptation from Marquardt (1963), elaborated in Bevington and Robinson (1992) and returns the χ 2 value, as goodness of fit parameter, a seven element vector containing, among the others, the widths of the Gaussians in the X and Y direction, the ellipse center location in pixel coordinate system, and, if present, the tilt angle of the ellipse in the clockwise direction from the X axis in radians. In Figure 4 (top) the best fitting surface (green) over-plotted to the surface representation of BO (yellow), imaged in Figure 3, is shown. Contour of the normalized values are plotted on the Z=0 plane. At the bottom of Figure 4 the fitting ellipses are overlaid to the rectilinearly gridded images of Figure 3. The U ellipse parameters actually refer to X (sample) and Y (line) couples in pixel coordinate system. Thus couples of longitude and latitude values can be associated to each geolocated pixel and corresponding am, bm path  has not been corrected. Sun illumination and viewing angles grow from right side to left side of the image: a) the structure at 373 nm appears as an elliptical shape encircled by a darker ring; another smaller oval is visible on the bottom-right corner of the image. b) at 752 nm the vortex shows many well-contrasted dynamic details; the previous smaller oval corresponds to the brightest feature at this wavelength. c) the deepest methane band at 890 nm enhances the high feature, located on the top of the vortex; the smaller oval appears darker (deeper) than its background. lengths can be calculated for the fitting ellipse in metric units (see the eccentricity equation ahead).
This operation implies a deformation of the original shape in pixel coordinates, connected to the warping that image undergoes once projected. The fitting ellipse axes have been converted in their corresponding path lengths through another code of the IDL library that returns the rhumb line -the line crossing all meridians of longitude at the same angle-that connects two points on a sphere, known the latitude and longitude of the two points and the radius of the sphere. As Saturn is not a spherical body, the local radius of curvature of the meridian passing for the vortex center had to be calculated in place of the sphere radius. The BO shape has been estimated through the eccentricity value in path lengths metric units, or:  where am and bm are the major and minor half path lengths in metric units. The size corresponds to the U area values. The shape tilt on the equator line is calculated as the regression coefficient of the ellipse major axis, or on the three couples of longitude (x), latitude (y) values of the U center and the two vertices of the largest pathlength.
In Figure 5 the U ellipses imaged in Figure  4 (bottom) in image plane coordinate system are visualized on the map of the observations in planetocentric coordinate system. Size, shape, areas and tilt are given in Table 2.
Errors on the BO size, shape and dimension are derived from the kilometric pixel size reported in Table 1. The kilometric pixel error is calculated as the diagonal half value of the square region around the pixel center. Errors connected to the fit parameters are negligible compared to the size error. The eccentricity (shape) error is obtained from the propagation of the path lengths uncertainties. Analogous computation has been carried out to evaluate the error on the area.
Latitude and longitude drift has been estimated following the BO center coordinate dis-   (Vasavada et al., 2006;Fletcher et al., 2012), some conclusions can be draw. Then we have interpolated the five points of the ellipse centers in latitude and longitude with polynomials fit, where the x axis is a Julian day's grid starting from the first observation of Table 1. The first derivatives of these curves have been also calculated on each grid point to obtain a propagation rate trend and finally the average on the whole time interval has been assumed as drift speed value. BO results to have a negligible mean propagation rate of 7.16 * 10 −3 ± 6.66 * 10 −3 deg/day toward the North. The mean propagation rate in westward direction is instead much faster being 0.75± 0.25 deg/day, corresponding to 8.9± 3.0 ms −1 in System III reference frame, in agreement with the results of Sayanagi et al. (2013). It is worth to observe that the only ISS-WAC observation used in this report is not decisive in the drift speed determination, however confirming the tendency of the BO to propagate in clockwise direction. In Figure 6 the BO center coordinates (red points), the best fit of the displacement (red line) and the first derivative of the best fit displacement (black dashed line), are reported both for the latitude (top frame) and for the longitude (bottom frame). The average propagation rates for the meridional (top) and zonal (bottom) components are depicted as black solid lines. Watching at the first derivative behavior, the BO meridional displacement would seem undergoes an acceleration with the time, but the low values associated (left y axis) induce to consider it stationary in latitude (top frame). Zonal acceleration of BO (bottom frame) shows instead a fast raising in the first 100 days to reach a maximum around the day 140, followed by a slow decrease. This is in line with a dynamic connected to the outburst interpretation, where the strong initial energy of the diffusive wake is damping with the time.
Summarizing the results reported in Table  2, the BO area seems to undergo a shrinking up to January 2012 (Figure 7). The shape evolution does not show a similar linear behavior, oscillating apparently in a random way during the first half of the year, as well as the tilt, and approaching a circular form at the end of the period. On the other hand the fit has been done on the shape of the image plane, thus a sensible uncertainty on the effective contour have to be considered as a bias.
In conclusion the mathematical procedure setup to analyze the life cycle of the top cloud over a vortex in the Saturn's giant storm on 2010 has been successfully checked. Dimensions are very different from those reported in Sayanagi et al. (2013), but this difference is probably due both to the different used spectral frames and to the considered cloud features. However size shrinking, direction and speed drift are confirmed. Future work is planned to apply this technique to other spectral frames with the aim to produce a 3D model of the vertical structure, which would help to constrain its dynamic origin.