PRACTISE – Photo Rectification And ClassificaTIon SoftwarE ( V . 1 . 0 )

Terrestrial photography is a cost-effective and easy-to-use method for measuring and monitoring spatially distributed land surface variables. It can be used to continuously investigate remote and often inaccessible terrain. We focus on the observation of snow cover patterns in high mountainous areas. The high temporal and spatial resolution of the photographs have various applications, for example validating spatially distributed snow hydrological models. However, the analysis of a photograph requires a preceding georectification of the digital camera image. To accelerate and simplify the analysis, we have developed the “Photo Rectification And ClassificaTIon SoftwarE” (PRACTISE) that is available as a Matlab code. The routine requires a digital camera image, the camera location and its orientation, as well as a digital elevation model (DEM) as input. If the viewing orientation and position of the camera are not precisely known, an optional optimisation routine using ground control points (GCPs) helps to identify the missing parameters. PRACTISE also calculates a viewshed using the DEM and the camera position. The visible DEM pixels are utilised to georeference the photograph which is subsequently classified. The resulting georeferenced and classified image can be directly compared to other georeferenced data and can be used within any geoinformation system. The Matlab routine was tested using observations of the north-eastern slope of the Schneefernerkopf, Zugspitze, Germany. The results obtained show that PRACTISE is a fast and user-friendly tool, able to derive the microscale variability of snow cover extent in high alpine terrain, but can also easily be adapted to other land surface applications.


Introduction
Oblique terrestrial photography has become a more and more frequently used observation method in various research disciplines, such as vegetation phenology (Richardson et al., 2007;Ahrends et al., 2008;Crimmins and Crimmins, 2008;Migliavacca et al., 2011), land cover studies (Clark and Hardegree, 2005;Zier and Baker, 2006;Roush et al., 2007;Michel et al., 2010) and volcanology (Major et al., 2009).Here, we focus on glaciology and snow hydrology where, for example, investigations of the snow albedo on glaciers were realised by Corripio (2004), Rivera et al. (2008) and Dumont et al. (2009).For a comprehensive overview of snow fall interception of vegetation, glacier velocity and snow cover mapping, we refer to Parajka et al. (2012).Terrestrial photography is used for these monitoring applications with an increasing frequency.This has to be attributed to the advancements in digital photography and in off-grid power supply.In addition, it is related to the fact that field campaigns and satellite-based remote sensing have limitations due to the prevailing weather conditions and the complexity of mountainous terrain (Klemes, 1990).Terrestrial photography offers an easy-to-use and inexpensive opportunity to monitor spatially distributed land surface characteristics, even in remote areas.
With the increasing availability of cost-effective highresolution digital cameras and high-resolution digital elevation models, new tools can be developed to observe and map the patterns of land surface variables such as the spatial distribution of the snow cover in mountainous terrain.The main   challenge for spatially distributed monitoring, however, is the georeferencing of a 2-D photograph to the 3-D reality.Tools, developed by Aschenwald et al. (2001) and Corripio (2004), addressed this problem by projecting the DEM to the camera image plane to establish a link between the photograph and the real world.Aschenwald et al. (2001) used a photogrammetric approach that needs various ground control points (GCPs) for the georectification process.This is, however, unfavourable in remote, mountainous terrain where the derivation of GCPs can be time-consuming and costly.Additionally, the integrated optimisation procedure in their approach only optimises the camera target position T , i.e. the centre position of the photograph, whereas all other parameters remain fixed.It should be noted here that T is known as the principal point in photogrammetry.The georectification method applied in Corripio (2004) is based on an animation and rendering technique by Watt and Watt (1992).This method only needs one GCP (T ), but 13 camera parameters have to be set.If these parameters are not accurately measured, they have to be manually corrected by changing them in an iterative way, which is unfavourable if extensive time series have to be processed.
The Photo Rectification And ClassificaTIon SoftwarE (PRACTISE) introduced here is based on the approach of Corripio (2004) but has been improved and extended by additional model features.We use slightly different formulations for the calculation of the 3-D rotation and projection.Even more importantly, several new optional routines are implemented in PRACTISE.This includes the dynamically dimensioned search (DDS) algorithm (Tolson and Shoemaker, 2007) to automatically identify the camera location and orientation using GCPs if the exterior and interior orientation parameters are not precisely known.Additionally, a view-shed algorithm (Wang et al., 2000) was integrated that simplifies and hastens the necessary visibility analysis by computing the viewshed directly without the additional step of using a geoinformation system, as is needed when using other georectification tools.PRACTISE also differs from existing software packages because it contains an automatic and a manual snow classification algorithm, and because a batch mode is implemented, i.e. several images can be classified in one program evaluation.As stated above, the routines described here are optional and the user selects the routines depending on the task and the available data.If, for example, the exact camera location and orientation of a photograph is known, the DDS optimisation with the need for additional GCPs can be omitted.By contrary, the DDS routine is absolutely necessary for the georectification procedure if the parameters are not precisely known.The strength of PRACTISE is that the new features form a flexible, fast and user-friendly processing tool for analysing spatially and temporally distributed land surface variables.A further strength is that the Matlab source code is freely available, and even though it is designed to classify the snow variability in mountainous terrain, it can be easily adapted to other fields of research, such as greenness indexes in phenology (Richardson et al., 2007;Ahrends et al., 2008;Crimmins and Crimmins, 2008;Migliavacca et al., 2011).

Data
The test area for PRACTISE is located near the Zugspitze mountain in the Alps (located in Bavaria, Germany, Fig. 1a).A common single lens reflex camera (SLR, Canon EOS 550D, Canon EF 17-40 mm f/4I USM objective lens) was installed at 2665 m a.s.l. at the Environmental Research Station Schneefernerhaus (UFS, 5th floor) which is located on the south slope of the Zugspitze.The camera was oriented towards the test area, on the northeast facing slope of the Schneefernerkopf summit (211 000 m 2 , Fig. 1b).The skiing area on the glacier was excluded.During daylight, hourly images were taken from 10 May 2011 to 2 March 2012 (307 days).Some technical problems with the automatic timer reduced the number of days with available photographs to 245 (2061 photographs).The hourly frequency, however, increased the probability that at least one suitable image would be obtained per day, which resulted in about 180 days with potentially suitable photographs unaffected by weather and lighting conditions.
PRACTISE requires as inputs a DEM raster and the exterior orientation parameters of the camera: the camera position C, the camera target position T and the roll φ of the camera.The latitude and longitude positions of C and T are sufficient as the altitude is taken from the corresponding DEM pixel during the computing process.If necessary, a camera offset o (installation height above the surface) is added to the altitude of C, the combined altitude being referred to as C o .Additionally, interior orientation parameters of the camera are necessary, such as the focal length f , as well as the sensor (CCD or CMOS) dimensions: height h and width w.We note here that lens distortions which can be significant are not taken in account in PRACTISE.Therefore, a high-quality objective lens was chosen that is known to have almost no distortion.
We selected three photographs to show the functionality of PRACTISE for different weather conditions and snow cover extents.The camera image taken on 11 May 2011 at 08:15 CEST (Fig. 1b) represented the start of the ablation period in spring under clear sky conditions.The photograph from 16 August 2011, 11:05 CEST, was recorded under clear sky conditions in summer with almost no snow in the investigation area, whereas the photograph from 17 February 2012 at 15:07 CET described cloudy conditions directly after a snowfall event in winter.
The input given for the georectification and the classification of the photographs is presented in Table 1.All camera dependent parameters were taken from the user manual of the Canon camera system.Using the best resolution (17.9 Mpx), the pixel dimensions of the photographs are vertically 3456 px (N v ) and horizontally 5184 px (N h ).The latitude and longitude positions of C and T were visually derived from an official orthophoto from September 2009 with a spatial resolution of 0.2 m provided by the Bavarian State Office for Survey and Geoinformation.The UFS building where the camera is located is clearly identifiable in the orthophoto, while the coordinates of T are estimated by comparing the orthophoto with a photograph from September 2011.Other techniques to obtain the coordinates of C and T might also be possible, for example with a standard GPS device.The parameters o, φ and f were estimated af- The available input data for the camera location and orientation is subject to considerable uncertainty as it was not accurately measured using for example a differential GPS system.Furthermore, the camera was moved between the images.The DDS optimisation utilising GCPs was applied to improve the exterior and interior orientation parameters of each photograph.The GCPs of each photograph were determined by using the orthophoto in combination with the DEM for the longitude, latitude and altitude as well as the photograph with the row and column information.

Model routines
PRACTISE is programmed in Matlab and divided into four modules that are presented in the following sections.The partitioning of the software in different routines provides a maximum of flexibility as the user can decide depending on the task and available data which features are necessary and have to be activated or if a new routine has to be implemented.In the default case, the camera location and orientation is precisely known.PRACTISE starts with the viewshed generation (Sect.3.1).Subsequently, the georectification procedure is applied (Sect.3.2) and finally the snow classification is executed (Sect.3.4).In our study, however, the exterior and interior orientation parameters of the camera are estimated.Hence, all routines are activated.In this case, PRACTISE begins by assessing the accuracy of the GCPs (Sect.3.3) where it utilises the georectification routine (Sect.3.2) to compute the deviations between the georeferenced and real positions of the GCPs.Next, the DDS algorithm reduces the positional inaccuracy by optimising the camera parameters (Sect.3.3).From this point, the default procedure is followed, as described above.In each section, we will show the processing steps based on the photograph of 11 May 2011 at 08:15 CEST (Fig. 1b).

Viewshed
In a first step, PRACTISE identifies the pixels of the DEM which are visible from the camera location.This is necessary because pixels of the digital image can only be attributed to those DEM pixels.Note that the spatial resolution of the DEM determines the detail of the results.The implemented viewshed calculation is an optional feature that can be bypassed if a viewshed is externally provided, for example from geoinformation software.
The viewshed generation is based on the "reference planes" concept (Wang et al., 2000) and requires a DEM raster and the camera position C o .By definition, only the horizontal centres of DEM pixels are utilised in the visibility analysis and the origin of the raster grid is in the northwestern (NW) corner.Indices i and m refer to the row positions of DEM pixels, and indices j and n indicate the column positions.
The viewshed calculation is divided into eight sectors based on the compass directions N, NE, E, SE, S, SW, W and NW.At first, the elevation of the DEM coordinate system is modified by setting the elevation of C o to zero.The normalised camera position simplifies the plane generation and is referred to as s i,j (Fig. 2a).The algorithm starts the visibility analysis at the DEM pixels in the second ring and proceeds stepwise to the cells of the outer rings.All pixels in the first ring are assumed to be visible, since no obstacles to s i,j are evident.The general functionality of the method is shown by using the example of the west-northwest (W-NW) sector (shaded area in Fig. 2a and b).
Three pixel values define the plane which builds the criteria of visibility (Z) for the destination point d m,n .These pixels are the normalised camera position s i,j as well as the neighbouring pixels r m,n+1 and r m+1,n+1 .Both, r m,n+1 and r m+1,n+1 , lie on the adjacent inner ring of d m,n , i.e. the third ring in Fig. 2b.Additionally, these two points have the shortest distance to s i,j and to d m,n on that ring.The values of r m,n+1 and r m+1,n+1 represent the maximum height of either the normalised elevation at this raster position or, in relative terms, higher obstacles in the already calculated inner rings in between to s i,j .
Z is then derived as follows:  The "reference planes" concept of Wang et al. (2000) is evaluated subsequently from the inner to the outer rings and is shown for an example in the 4th ring of the greyshaded W-NW sector.The normalised camera position s i,j , as well as the neighbouring pixels r m,n+1 and r m+1,n+1 (third ring) create a plane that checks if the pixel with the normalised elevation value d m,n is visible.In this case, d m,n is visible as the plane height Z at row m and column n is lower (adapted from Wang et al., 2000).
The calculation for the main directions is simplified since the reference plane (Eq. 1) can be reduced to a "reference line" (Eq.2).This is shown for the NW diagonal: A pixel is considered as visible if In this case the value of d m,n is assigned to r m,n for further calculations in the adjacent outer ring, otherwise the pixel is invisible and r m,n is set to the value of Z.The next visibility check will be executed at d m−1,n (Fig. 2b).Other directions and sectors are calculated in a similar way.The algorithm of Wang et al. (2000) was developed to generate a 360 degree viewshed.Assuming a central projection of the camera lens, we use the viewing direction as well as the horizontal and vertical field of view and thus only compute the areas depicted on the photographs.Here, we additionally need the camera target position T and the interior orientation parameters of the camera f , h and w.The viewing direction is set by connecting C o and T .The interior orientation parameters are necessary to calculate the corresponding horizontal and vertical field of view.A maximum vertical viewing angle α v to the viewing direction can be calculated as follows: The maximum horizontal viewing angle α h of the photograph is calculated by replacing the height h by the width w in Eq. ( 3).The vertical or horizontal orientation of a camera image might be different to the real world vertical or horizontal orientation due to φ.
Figure 3 shows the viewshed in this case study.

Georectification
PRACTISE uses an animation and rendering technique to georectify the visible DEM pixels (Watt and Watt, 1992).21 Fig. 4. The principle of the georectification procedure is as follows: at first, the mountain massif in the real world coordinate system (XYZ W , black) is translated and rotated to the camera coordinate system (XYZ C , blue).Then, the 3-D mountain landscape is projected to a 2-D virtual camera image utilising the central projection of the camera lens (adapted from Corripio et al., 2004).
The principle behind the georectification process is illustrated in Fig. 4. The camera produces a 2-D representation of the 3-D landscape.The oblique and two-dimensional image lacks depth information: therefore a direct back-calculation of the 2-D information into a 3-D landscape is impossible.
A method for calculating it is to generate a 2-D virtual camera image of the DEM while conserving the real world position of any pixel.The RGB (red, green, blue) values of the camera can then be assigned to the virtual 2-D image.
Afterwards, any pixel with the attached RGB information is retransformed to its real world position.
The georectification is shown for a single DEM pixel, whereas all visible pixels are successively processed in the same way.Given the fact that the pixel is visible from C o , its centre point coordinates P w are derived and saved in a vector: The transformation of the real world coordinates of the DEM into the camera coordinate system is achieved by a translation of the origin of the coordinate system to the camera position C o and a subsequent multiplication of the translated pixel coordinates with a rotation matrix:  22 Fig. 5.The mathematical components of the translation and rotation of the real world coordinate system (XYZ W , black) can be derived using vector calculus.The translated real world coordinate system (XYZ T , red) is determined by setting C o as coordinate system origin.The connection line from C o to T forms the vector of the viewing direction which is subsequently normalised (N ).The unit vector U is derived by the cross product of N and the unit vector of N xy (green) where N xy is the projection of N to the XY T plane.The directions of the camera coordinate system (blue) are spanned by N , U and V where V is the cross product of U and N (adapted from Corripio, 2004).
The unit vectors U , V and N describe the axis of the new camera coordinate system (Fig. 5), where N points in the viewing direction.Vectors U and V are the horizontal and vertical axis of the camera system and create a plane that is parallel to the image plane (Figs. 4, 5).The calculation of N is performed on the basis of the real world coordinates of C o and T : Following Corripio (2004), we use cross products to calculate U and V (Fig. 5), if N z = 0: where N xy = [N x , N y , 0].We extend the calculation to the situation where N = N xy , i.e.N z = 0.In this particular case, we set V = [0, 0, 1], and calculate U by computing the cross product of V and N.
In the case of the camera not being completely levelled, an additional rotation of the coordinates around N is required: where the roll φ is defined from 0 • to ± 90 • , where the positive values turn the U -V plane in the viewing direction clockwise, and the negative values turn it anticlockwise.The last step of the georectification is the projection of the rotated coordinates P cr to the image plane p.The three coordinate values of the DEM pixel determine the position in the camera space, where P crx as well as P cry hold the horizontal and vertical information and P crz the depth information.In contrast to Corripio (2004), we reduce the 3-D problem (P cr to P p ) to two 2-D problems: a horizontal (P crx to P px ) and a vertical (P cry to P py ) one.We solve each of them in two steps.At first, we directly apply the intercept theorem to calculate the horizontal (and vertical) component of the photograph at the CCD sensor plane s: where f is the focal length.P sy is calculated by replacing P crx with P cry .It should be noted that the intercept theorem is applicable here, as T is located at [0, 0, P crz ] in the camera coordinate system and thus lies in the centre of the photograph.As a second step, P sx and P sy are scaled to the image plane p using the number of pixels N h and N v of the photograph in the horizontal and vertical directions: where w is the camera sensor width.P py is computed in the same way but under usage of N v and the CCD height h.Both, the photograph and the projected DEM are now in-plane.The last step is to shift the origin of the virtual camera image from to the origin of the photograph [0, 0].This is necessary as the photograph origin lies at the upper left corner while the projected DEM coordinate system is centred in the photograph.The overlay of the images facilitates the direct extraction of the RGB values for the classification which can be directly retransformed to the raster format of the DEM.
Figure 6 shows the overlay of the georectified DEM pixels and the photograph.

GCP accuracy assessment and DDS optimisation
PRACTISE offers an optional feature to enhance the exterior and interior orientation parameters of the camera used in the georectification procedure if the camera parameters are not Fig. 6.The georectification of the visible DEM pixels (Fig. 3) is superimposed with cyan dots on the corresponding photograph.precisely known.In that case, GCPs are required to determine and to reduce the positional inaccuracy of the virtual camera image to the photograph.The root mean square error (RMSE) is used as an error metric.
We implemented a global optimisation approach, the dynamically dimensioned search (DDS) algorithm (Tolson and Shoemaker, 2007), to minimise the displacement between the georeferenced and real locations of the GCPs.We selected this technique because Tolson and Shoemaker (2007) state that, for calibration problems between 6 and 30 dimensions and with a limited number of function evaluations (1000 to 10 000), it produces equally good or even better results than the frequently used shuffled complex evolution (SCE) optimisation.The general procedure of the implemented DDS algorithm is shown in Table 2.
Within the optimisation procedure, seven decision variables are optimised: the latitude and longitude of C and T , the camera offset o, the roll φ, and the focal length f .The inclusion of f in the optimisation is necessary as the actual and nominal focal length of a camera lens will probably differ slightly.The initial estimates of the decision variables x 0 are taken from the original input (Table 1).Additionally, the user has to define the upper and lower boundaries (x max and x min ) that span the range of reasonable values (Table 3).Finally, the maximum number of function evaluations m is specified.Six GCPs are used in the DDS optimisation example.
The algorithm starts with the georectification of the GCPs using x 0 and creates an initial x best = x 0 .Then, x new is randomly generated (Table 2) and if the recalculation of the RMSE(x new ) results in a lower RMSE than RMSE(x best ), x best is updated with x new .The optimisation procedure stops when the number of iterations is equal to m and subsequently the georectification of the DEM with the best camera orientation (Table 2) starts.In this example, m = 3000 as no large improvements have been observed with more iterations.At least one recalculation is recommended to verify that the global optimum was found.
Figure 7 depicts the correct position of the six GCPs (green crosses) in comparison to the georectification of the GCPs before and after the DDS optimisation (red circles and dots).

Classification
Here we focus on the classification of snow cover even though the investigation of other land surface variables is possible and only needs slight adaptions of the respective routine.Two classification routines of different complexity can be used.The first is based on threshold values, which have to be manually derived by analysing the RGB values of the snow cover and of the surrounding environment in the photograph.The second is an automatic snow cover classification routine (Salvatori et al., 2011), that has not been used in quantitative snow cover mapping before but is able to hasten the classification, in particular of long time series.
The manual classification assigns snow to pixels with RGB values above certain thresholds.The threshold values can be in between 0 and 255 if 8 bit data is used and are often around 150 for all bands.The predefined thresholds vary from image to image as the lighting conditions change and as we want to classify fresh snow (pure white), as well as old snow which turns grey with time.Snow is generally approximately equally reflective within the RGB bands, while the reflectance values of, for example, light-coloured bare rocks are significantly lower in the blue band.Hence, we introduced a test that verifies if the spread between the RGB values of one pixel does not exceed a specified threshold (for example 10).
The automatic classification of Salvatori et al. (2011) incorporates a statistical analysis of the image by using a DN (digital number) frequency histogram (Fig. 8a).The algorithm uses the blue band exclusively because of the assumption that it is representative for the other bands with respect to snow.In the presence of snow, the histogram usually shows a bimodal distribution.The first local minimum over or equal to 127 is selected as the snow threshold (Fig. 8a).The DN frequency histogram has to be smoothed for this analysis by using a moving average window of 5.This is done for removing single outliers, which might be mistakenly interpreted as local minima.Salvatori et al. (2011) defined the size of the moving window as well as of the minimum histogram threshold of the blue band, on the basis of about 300 images.The resulting classification is shown in Fig. 8b.
The structure of PRACTISE also allows for an inclusion of already classified images and of other routines.For example, Hinkler et al. (2002) present a calibrated index similar to the normalised-difference snow index in satellite remote sensing (Dozier, 1989) to identify snow cover and areas free of snow, while Schmidt (2007) uses manually determined thresholds and additional masks of shadows, vegetation and topographic features.In the studies of Corripio (2004) and Corripio et al. (2004), the albedo of glacier and snow surfaces is calculated using an atmospheric transmittance model.Algorithms for the investigation of other land surface variables are likewise possible.The implementation of any existing or self-programmed routines in PRACTISE can also be accomplished with limited programming skills.

Results and discussion
The functionality of PRACTISE will be demonstrated using the test area of Schneefernerkopf, Zugspitze on the basis of three photographs which are hereinafter referred to as the May (11 May 2011 at 08:15 CEST), the August (16 August 2011, 11:05 CEST) and the February (17 February 2012 at 15:07 CET) images.All of the described routines are used to compute the snow cover extent.The DDS optimisation is necessary as the exact camera location and orientation was not measured but estimated from an orthophoto.Further, the camera was slightly moved between each photograph.The runtime per photograph is about 40 s using computing power similar to an Intel Pentium 4 with 3 GHz.
In the DDS optimisation, we found an initial RMSE of 67.82 px between the GCPs and the control points within the May photograph (Table 3, x 0 ).The error could be reduced to 4.42 px by using the optimised input (Table 3, x best ).The RMSE after the optimisation procedure for the August and February images are 6 and 5.49 px, whereas the initial error values are 43.45 and 92.91 px, respectively.The comparison of the RMSE values illustrates that the positional accuracy of the optimised input is at least seven times higher than the initial input.The mean RMSE of these three photographs (5.30 px) corresponds to 0.79 m for the mean distance of 1044.46 m between the camera position and the GCPs and is thus smaller than the spatial resolution of the DEM (1 m).
The visual investigation of the automatically classified photographs showed qualitatively a good agreement between automatically classified and visually observed snow covered areas (Fig. 9a-c, enlarged views of the investigation area).The high quality of the classification applies to both clear sky conditions in the May image (Fig. 9a) and cloudy conditions in the February image (Fig. 9b).In the August photograph (Fig. 9c), limitations in the classification with respect to light-coloured bare rock could be observed.A small test area (11 701 m 2 , black box in Fig. 9c) was selected within  Tolson and Shoemaker, 2007).
Step Step 4 FOR j = 1, ..., J decision variables in {N }, perturb x best j using a standard normal random variable N (0, 1), reflecting at decision variable bounds if necessary: x new j = x best j + σ j N (0, 1), where σ j = r(x max j − x min j ) IF x new j < x min j , reflect perturbation: Step 5 Evaluate the investigation area where visually no snow could be detected.However, the automatic classification routine mistakenly classifies 477 m 2 of limestone as snow, which corresponded to a relative error of 4.1 %.The August photograph was also processed using the manual classification routine (Fig. 9d).The thresholds of the RGB bands were identical to the automatically derived Table 3. Vectors of the DDS optimisation example with 3000 iterations: x 0 , x max , x min and x best .The values are in m except noted otherwise.x max of C is set to the values of x 0 as the UFS building is represented in the DEM by a plateau.Hence, we confine the optimisation directions to stay at the edge or in front of the building.The latter needs a large camera offset to obtain the height of the 5th floor of the UFS.classification threshold (169) in Fig. 9c.The maximum allowed spread between the three RGB values of one pixel was 10.Qualitatively, the visual investigation of the overlay of the photograph and the classification shows a good match for the investigation area.We investigated again the same small test area (black box in Fig. 9d).The misclassification is reduced to 100 m 2 (0.9 %) in comparison to the automatic classification, due to the light-coloured bare rock reflecting the blue band significantly more weakly than the red and green bands.
The resulting snow cover maps of the three photographs are depicted in the Fig. 10a-c using the classification of Fig. 9a, b and d, respectively.We compared the derived snow cover extent in Fig. 10a, to the DEM, respective to the slope in more detail.More than 90 % of the areas free of snow on this date are located in steep terrain with slope angles above 35 • (without figure).With the last snowfall being on 3 May 2011, this is reasonable due to gravitational snow redistribution (Bernhardt and Schulz, 2010).The snow cover extents in the investigation area (black dotted line) Fig. 10.The maps depict the resulting snow cover extent of the Fig. 9a (a), 9b (b), and 9d (c).The black dashed line outlines the investigation area at the Schneefernerkopf.We want to note here that the small test area (black box in Fig. 9d) is not shown in (c).
26 Fig. 10.The maps depict the resulting snow cover extent of the Fig. 9a (a), 9b (b), and 9d (c).The black dashed line outlines the investigation area at the Schneefernerkopf.We want to note here that the small test area (black box in Fig. 9d) is not shown in (c).are in accordance with the time of the year, and amount to 94 000 m 2 on 11 May 2011 (Fig. 10a), 122 000 m 2 on 17 February 2012 (Fig. 10b) and 13 000 m 2 on 16 August 2011 (Fig. 10c).
The present results reveal that PRACTISE, with its different features and its flexibility, is an efficient software tool to produce temporal and spatial high-resolution snow cover maps.All methods used are well-established and the optional routines can be selected by the user depending on the available data and the task.We have shown here that the DDS optimisation as well as the classification routines produce high quality results for the three investigated photographs.The accuracy assessments of all three images are better than the spatial resolution of the DEM.Thus, the DDS optimisation of the interior and exterior orientation parameters makes the software very valuable for the analysis of photographs where the camera parameters are only imprecisely known, for example for extensive time series where camera movements are a problem.The combined use of DDS optimisation and viewshed routine additionally hastens the georectification procedure in our study as with each new camera location, a new viewshed is needed.The automatic snow classification of Salvatori et al. (2011) works well in most cases without the need for calibration or the manual determination of thresholds for different weather situations and snow cover patterns.The manual classification routine provided in PRACTISE can be used as an alternative under unfavourable conditions.Although the automatic classification represents a promising approach, the presented classification results also confirm the well-known limitations in the snow classification using the visible spectrum (0.4-0.7 µm).Shadows are another possible source of uncertainty; however, they do not have a great effect on the snow cover mapping here as the recording time was able to be controlled and was adjusted to a minimum of shading by choosing the day time (Dozier, 1989;Winther and Hall, 1999;Schmidt, 2007;Salvatori et al., 2011).
The fast and easy processing capabilities of PRACTISE might help to increase the efficiency of terrestrial photography either in validating spatially distributed snow hydrological models (Lehning et al., 2006) or in statistically analysing snow patterns influenced by the topography (Lehning et al., 2011).Future studies using PRACTISE will test the comparability of SLR images to other methods of snow cover detection and include long-term studies.A further topic of research will be the development of an automatic classification algorithm that is less prone to misclassifications of snow in digital camera images caused by clouds, shadows or lightcoloured bare rock.The versatility and the opportunity to comfortably georeference especially large time series of photographs in PRACTISE makes the software also attractive for other research disciplines.In particular, the aforementioned example of phenological greenness indexes, as well as the observation of land surface temperatures using thermal infrared cameras might be interesting fields of application.

Fig. 1 .
Fig. 1.(a) The test site of PRACTISE is located at the Schneefernerkopf which is situated in southern Germany, at the border to Austria (right frame).The DEM depicts the camera location and the field of view of the camera.(b) The installed digital camera system records hourly photographs of the investigation area, the north-eastern slope of the Schneefernerkopf summit (upper central area). 18

Fig. 1 .
Fig. 1.(a) The test site of PRACTISE is located at the Schneefernerkopf which is situated in southern Germany, at the border to Austria (right frame).The DEM depicts the camera location and the field of view of the camera.(b) The installed digital camera system records hourly photographs of the investigation area, the north-eastern slope of the Schneefernerkopf summit (upper central area).

Fig. 2 .
Fig. 2. (a) The viewshed calculation is divided into eight sectors based on the compass SE, S, SW, W and NW (black lines) from the point s i,j .The DEM pixels are attributed to 3, 4, ..., black dotted lines) depending on the pixel distance to s i,j .(b) The "reference plan et al. (2000) is evaluated subsequently from the inner to the outer rings and is shown for a ring of the grey-shaded W-NW sector.The normalised camera position s i,j , as well as the r m,n+1 and r m+1,n+1 (third ring) create a plane that checks if the pixel with the norma d m,n is visible.In this case, d m,n is visible as the plane height Z at row m and column fromWang et al., 2000).

Fig. 2 .
Fig. 2. (a)The viewshed calculation is divided into eight sectors based on the compass directions N, NE, E, SE, S, SW, W and NW (black lines) from the point s i,j .The DEM pixels are attributed to a certain ring (1, 2, 3, 4, ..., black dotted lines) depending on the pixel distance to s i,j .(b) The "reference planes" concept ofWang et al. (2000) is evaluated subsequently from the inner to the outer rings and is shown for an example in the 4th ring of the greyshaded W-NW sector.The normalised camera position s i,j , as well as the neighbouring pixels r m,n+1 and r m+1,n+1 (third ring) create a plane that checks if the pixel with the normalised elevation value d m,n is visible.In this case, d m,n is visible as the plane height Z at row m and column n is lower (adapted fromWang et al., 2000).

Fig. 3 .Fig. 3 .
Fig. 3.The optional viewshed feature of PRACTISE computes the visible pixels (cyan) using the corresponding camera location and orientation.

Fig. 4 .
Fig.4.The principle of the georectification procedure is as follows: at first, the mountain massif in the real world coordinate system (XYZW, black) is translated and rotated to the camera coordinate system (XYZC, blue).Then, the 3-D mountain landscape is projected to a 2-D virtual camera image utilising the central projection of the camera lens (adapted fromCorripio et al., 2004).

Fig. 5 .
Fig.5.The mathematical components of the translation and rotation of the real world coordinate system (XYZW, black) can be derived using vector calculus.The translated real world coordinate system (XYZT, red) is determined by setting Co as coordinate system origin.The connection line from Co to T forms the vector of the viewing direction which is subsequently normalised (N ).The unit vector U is derived by the cross product of N and the unit vector of N xy (green) where N xy is the projection of N to the XYT plane.The directions of the camera coordinate system (blue) are spanned by N , U and V where V is the cross product of U and N (adapted fromCorripio, 2004).

Fig. 7 .Fig. 6 .
Fig. 7.The correct GCP positions are depicted as green crosses in the enlarged view of the photograph.The georectification using x 0 is shown by the red circles while the red dots illustrate the georectification after the DDS optimisation (m = 3000) using x best .

Fig. 7 .Fig. 7 .
Fig. 7.The correct GCP positions are depicted as green crosses in the enlarged view of the photograph.The georectification using x 0 is shown by the red circles while the red dots illustrate the georectification after the DDS optimisation (m = 3000) using x best .
solution x 0 = [x 1 , ..., x 7 ] Vectors of upper, x max , and lower, x min , boundaries for the 7 decision variables Maximum number of function evaluations m Neighbourhood perturbation size parameter r (0.2 is default) Step 2 Set counter to 1, i = 1, and evaluate RMSE at initial solution RMSE(x 0 ): RMSE best = RMSE(x 0 ), and x best=x 0 Step 3 Randomly select J of the D decision variables for inclusion in neighbourhood {N}: Calculate probability each decision variable is included in {N } as a function of the current iteration count: P (i) = 1 − ln(i)/ ln(m) FOR d = 1, ..., D decision variables, add d to {N } with probability P IF {N} empty, select one random d for {N}

Fig. 8 .
Fig. 8. (a) The automatic snow classification in PRACTISE creates a DN frequency histogram of the blue band values (blue-green bars) of the superimposed DEM pixel positions (Fig. 6).The distribution is smoothed with a moving average window size of 5 (black line) and the snow threshold (green line) is selected for the first local minimum beyond a DN ≥ 127 (red line).(b) In the overlay, all DEM pixels with a DN in the blue band in the range from the snow threshold to 255 are classified as snow (red dots), while all other pixels are assigned as no snow (blue dots).

Fig. 9 .
Fig. 9.The superimposition of the DEM pixels (red dots = snow, blue dots = no snow) over the corresponding and enlarged photograph are shown on the left for the automatically classified images under clear sky conditions in spring on 11 May 2011 at 08:15 a.m.(a), under cloudy conditions in winter on 17 February 2012 at 03:07 p.m. (b), and under clear sky conditions in summer on 16 August 2011, 11:05 a.m.(c), as well as for the manually reprocessed classification of the August image (d).On the right, the corresponding snow thresholds (green lines) are illustrated: 153 (a), 134 (b) and 169 (c).The manual snow classification threshold is 169 for all three RGB bands and 10 for the maximum-minimum test (d).The black box in (c) and (d) depicts a small test area in the investigation area at the Schneefernerkopf where where visually no snow could be detected although several pixels are classified as snow.

Fig. 9 .
Fig. 9.The superimposition of the DEM pixels (red dots = snow, blue dots = no snow) over the corresponding and enlarged photograph are shown on the left for the automatically classified images under clear sky conditions in spring on 11 May 2011 at 08:15 CEST (a), under cloudy conditions in winter on 17 February 2012 at 15:07 CET (b), and under clear sky conditions in summer on 16 August 2011, 11:05 CEST (c), as well as for the manually reprocessed classification of the August image (d).On the right, the corresponding snow thresholds (green lines) are illustrated: 153 (a), 134 (b) and 169 (c).The manual snow classification threshold is 169 for all three RGB bands and 10 for the maximum-minimum test (d).The black box in (c) and (d) depicts a small test area in the investigation area at the Schneefernerkopf where visually no snow could be detected although several pixels are classified as snow.

Table 1 .
Initial input data of PRACTISE for the test site Schneefernerkopf.The coordinates are referenced to the European Terrestrial Reference System 1989 (ETRS89) and UTM Zone 32T.The values are in m except noted otherwise.

Table 2 .
Working steps of the implemented DDS algorithm in PRACTISE (adapted from RMSE(x new ) and update current best solution if necessary: IF RMSE(x new ) ≤ RMSE best , update new best solution: RMSE best = RMSE(x new ) and x best = x new