Detection of ground fog in mountainous areas from MODIS ( Collection 051 ) daytime data using a statistical approach

The mountain cloud forest of Taiwan can be delimited from other forest types using a map of the ground fog frequency. In order to create such a frequency map from remotely sensed data, an algorithm able to detect ground fog is necessary. Common techniques for ground fog detection based on weather satellite data cannot be applied to fog occurrences in Taiwan as they rely on several assumptions regarding cloud properties. Therefore a new statistical method for the detection of ground fog in mountainous terrain from MODIS Collection 051 data is presented. Due to the sharpening of input data using MODIS bands 1 and 2, the method provides fog masks in a resolution of 250 m per pixel. The new technique is based on negative correlations between optical thickness and terrain height that can be observed if a cloud that is relatively plane-parallel is truncated by the terrain. A validation of the new technique using camera data has shown that the quality of fog detection is comparable to that of another modern fog detection scheme developed and validated for the temperate zones. The method is particularly applicable to optically thinner water clouds. Beyond a cloud optical thickness of ≈ 40, classification errors significantly increase.


Introduction
Cloud forests are tropical and subtropical forest ecosystems characterized by the frequent occurrence of ground fog conditions (Bruijnzeel et al., 2010).As they intercept water from cloud droplets and, due to their mostly wet canopy, have a decreased rate of transpiration, they play an important role as an ecosystem service provider increasing local water supplies (Mildenberger et al., 2009).Furthermore they are biodiversity hot spots with a high number of endemic species (Postel et al., 2005).This also holds true for the cloud forest of Taiwan (Hsieh, 2002).While cloud forest areas of Taiwan are already the subject of intensive research (cf., e.g., Mildenberger et al., 2009;Chu et al., 2012), Taiwan's cloud forest has never been completely mapped on a country-wide scale.The most comprehensive information about the extent of cloud forest in Taiwan available today is given by Li et al. (2013), based on the National Vegetation Database of Taiwan.Since these data are based on field surveys, they are highly reliable but only cover the area of 9822 plots (each with an area of 400-2000 m 2 ) distributed over the whole country.Due to the inaccessibility of Taiwan's mountainous areas, those plots are mainly located close to roads.
For a spatial-explicit mapping, the usage of remotesensing data seems advisable.As the occurrence of cloud forest depends on the heavy influence of ground fog conditions, it has been shown by Mulligan and Burke (2006) that it can be discriminated from other forest types by the application of a threshold on maps of the ground fog frequency.The potential of this approach for the application in Taiwan has been shown by Thies et al. (2015) using lowstratus frequency maps derived from Moderate-resolution Imaging Spectroradiometer (MODIS) data.Because no distinction has been made in this study between low-stratus clouds (including clouds without ground contact) and ground fog, the explanatory power of the presented low-stratus frequency maps is, however, limited.To obtain more significant ground frequency maps, an algorithm able to detect ground fog (defined here as any cloud with ground contact) in Taiwan from satellite data is necessary.Common techniques for ground fog detection are analyzed with respect to their applicability for Taiwan in Sect.2.1.A new method which is more suitable for mountainous areas is presented in Sect.2.2.

Existing approaches
In order to detect ground fog from space, commonly a plane cloud base is assumed.The height of the cloud base is compared to a digital elevation model (DEM).If it is equal to or below the terrain height taken from the DEM, ground fog conditions can be assumed.Mulligan and Burke obtain the height of the cloud base surface by modeling the lifting condensation level of clouds (the area of which is taken from High-resolution Infrared Radiation Sounder (HIRS) data) from monthly WorldClim data that are interpolated from station measurements (Hijmans et al., 2005).For a higher temporal resolution and a more precise cloud base, usually the height of the cloud top and the cloud geometrical thickness are retrieved from satellite data using the abstraction of a plane-parallel cloud geometry (not to be confused with plane-parallel cloud models, which are used in radiative transfer calculations).This is possible as ground fog conditions are commonly caused by stratiform clouds.The height of the cloud base can then be calculated by subtracting the cloud geometrical thickness from the cloud top height (cf.Fig. 1).
Different methods for the cloud top height retrieval, e.g., CO 2 slicing (Menzel et al., 1983), DEM extraction (Bendix and Bachmann, 1993), or a method recently presented by Yi et al. (2015), do exist.As they are not suited for low clouds or only work under certain conditions, they cause different problems when applied to ground fog detection over mountainous areas (cf.Table 1).Cloud top heights may also be calculated from the cloud top temperature retrieved from infrared window channels.Cermak and Bendix (2008) compare it to the temperature of surrounding land pixels.Under the assumption of a certain (negative) temperature lapse rate with altitude the height of the cloud top above ground can then be calculated.Besides the problem that the assumed temperature lapse rate may not be given in many cases, this approach is neglecting material properties such as albedo or heat capacity and therefore the different reactions of cloud and ground pixels to illumination or temperature changes.It has, however, been successfully tested in the temperate zones of the Earth (Cermak and Bendix, 2011) so that these problems do not generally seem crucial.In Taiwan the situation seems to be different, as shown in a MODIS scene containing ground fog in Taiwan (Fig. 2).The fog is limited by the surrounding terrain and its surface is relatively flat.Therefore its surface height should be at the height of the lowest visible land surface.According to the temperature lapse rate approach it should also have the same temperature.The latter is obviously not the case.Instead the temperature difference between the fog and the land surface is about 10 K, which would result in a height difference of at least several hundreds of meters if a reasonable lapse rate is assumed.
Different methods are also available for cloud geometrical thickness derivation (cf.Table 2).Simple approaches that use empirically derived relationships between the thickness and the liquid water path (LWP) of a cloud (Hutchinson, 2002) or its optical thickness (Minnis et al., 1997) lack precision because of oversimplification.More sophisticated approaches applied in ground fog detection schemes as pseudosounding (cf.eg.Chang and Li, 2002;Bendix et al., 2005) or an approach based on iteratively simulated LWPs used by Cermak and Bendix (2011) use more complex cloud parameterizations but rely on many assumptions regarding cloud microphysical properties and their vertical distribution within the cloud.These assumptions may be valid for radiation fog.The typical Taiwanese mountain fog, however, seems to be -based on field experience and information given by Li  The CO 2 absorption is very high for low levels (< 3 km) of the atmosphere.This results in bad signal-to-noise ratios for low level clouds such as fog.Bendix and Bachmann (1993) DEM extraction: for a fog entity that is horizontally restricted by the terrain the height of its outermost pixels can be read from a DEM.
Only possible for fog that is restricted by the terrain in its horizontal extent.

Yi et al. (2015)
The top height of radiation fog that is restricted in its vertical extent by an temperature inversion is equal to the base height of that inversion.
The base height of the inversion can be calculated from satellite data.

Highly experimental.
Works only for radiation fog.The fog occurrence in Taiwan, however, is mainly caused by moist air masses being uplifted by the Taiwanese mountains (Li et al., 2015).
For example, Cermak and Bendix (2008); Platnick et al. (2003) Under the assumption of a fixed negative temperature lapse rate or an atmospheric profile, the cloud top height can be calculated from the cloud top temperature.
Neglects material parameters.This is crucial for ground fog in Taiwan.
An assumed temperature lapse rate/profile might be very wrong for many scenes.As the liquid water path (LWP) is the column integration of the the liquid water content (LWC), the cloud thickness can be calculated from the satellite-retrieved LWP under the assumption of a fixed LWC for certain cloud types.
The LWC is vertically not constant.For thin clouds, however, a vertically constant LWC can be approximated.Minnis et al. (1997) The cloud thickness is calculated from the satellite-retrieved cloud optical thickness using empirical formulas.For clouds in different heights, different formulas are used.
Due to oversimplification, the approach can only be seen as a crude approximation.
For example, Chang and Li (2002); Bendix et al. (2005) Pseudosounding: measured albedos in different channels of the solar spectrum are compared to theoretical albedos that were simulated for clouds with different thicknesses using radiative transfer calculations and stored in lookup tables.
The thickness of the simulated cloud with the smallest deviation between its albedos and the measured albedos is then assumed for the real cloud.
For the radiative transfer calculations several assumptions about the cloud microphysics are necessary.These assumptions may be not true for Taiwanese fog clouds.Cermak and Bendix (2011) Clouds with different cloud thicknesses are iteratively simulated using a three-layer cloud model.The LWC of the simulated cloud is integrated over the cloud thickness in order to obtain the LWP.This theoretical LWP is compared to a satellite-retrieved LWP.If they match, the thickness of the simulated cloud is assumed for the real cloud.
Several assumptions about the cloud microphysics are necessary.These assumptions may be not true for Taiwanese fog clouds.
et al. ( 2015) -more of advective nature.While it has not been tested whether existing cloud thickness retrievals can be applied despite some inaccuracies, the main problem of an accurate cloud top height derivation would remain.

The new approach -theoretical preliminary considerations
The main problem in ground fog detection in Taiwan is the derivation of cloud top heights.A method that does not take this intermediate step and detects the cloud base directly instead seems to be suited to overcome it.
For an assumed plane-parallel cloud as shown in Fig. 1, the geometrical thickness is by definition the same for all parts of the cloud that do not touch the terrain.If that cloud is, however, cut off by the terrain in some parts (causing ground fog), its thickness would be reduced in this area.If it is additionally assumed that the cloud is horizontally homogeneous in its optical properties, the optical thickness of that cloud would correlate with its geometrical thickness.Therefore the optical thickness should be constant in the parts without ground contact.In the parts with ground fog it would be reduced (decreasing with an increasing height of the terrain).For such an assumed cloud, a simple filling algorithm could be applied on an image of its optical thickness in order to flag all the pixels with the same values (cf.Fig. 3).In the flagged pixels the cloud would not have ground contact (red, right side).The unmarked cloud pixels (shades of blue, right side) would have ground contact causing ground fog.
A real-world stratiform cloud is neither perfectly planeparallel nor horizontally homogeneous in its optical properties.It can be described by this model only to some degree.Therefore a simple filling approach is obviously not suitable for fog detection.The pixels without ground contact would still have different optical thickness values.It can, however, be assumed that no statistical relationship between the height of the terrain and the optical thickness is given in the parts of the cloud that are not touching the ground.In contrast, there is a negative correlation between the terrain height and the optical thickness for the pixels with ground fog as the geometrical thickness decreases with increasing terrain height.Without knowledge about the vertical distribution of a cloud's extinction coefficient a linear correlation cannot be assumed.It is only known that the optical thickness decreases with increasing DEM height.Such a trend can be detected using Spearman's rank correlation coefficient (ρ).
Let ρ below h be Spearman's rank correlation coefficient calculated from the DEM height and the remote-sensed optical thickness of all pixels of a cloud entity with a DEM height below a certain height h.If we assume perfect conditions (a perfectly plane-parallel cloud without any horizontal inhomogeneities and a perfectly retrieved -no sensor noise, for example -optical thickness), ρ below h should be 0 for all h at the level of or below the cloud base height.For all h above the cloud base height it should be below 0 and decrease with increasing h.Let ρ above h be the correlation between height and optical thickness calculated from all cloud pixels with a height above or equal to h.Under perfect conditions ρ above h should be −1 for all h at the level of or above the cloud base height.For all h below the cloud base height it should be above −1 and increase with decreasing h.Based on this, the height in which the difference ρ diff h = ρ below h − ρ above h is maximal can be defined as the base height of a cloud.This results in a small ρ diff p .(b) At the cloud base height, ρ above p is calculated from pixels with ground fog.It is therefore clearly negative.ρ below p is calculated from pixels without ground fog.Its value is therefore close to 0. This results in a high ρ diff p .(c) Below a cloud, ρ above p and ρ below p are calculated from pixels without ground fog.Therefore both values are near to 0. This results in a small ρ diff p .
As Fig. 4 shows, this theoretically assumed relation between optical thickness and terrain height can be observed for real clouds with ground contact.Above the measured cloud base height the optical thickness decreases with the height of the DEM as the cloud is cut off by the terrain here.Below the cloud base height, no relationship between optical thickness and DEM can be observed.In this real-world example, ρ above cloud base does, of course, not reach values as extreme as −1 and also ρ below cloud base is not exactly 0.
The domain of the scatter plot in Fig. 4 is relatively small on purpose.Further in the north -where the DEM values are decreasing -the optical thickness also decreases.Thus, the assumption of plane-parallelism and/or horizontal homogeneity is clearly not valid for the depicted cloud.On a more local level, however, the variability of cloud top height, cloud bottom height, as well as the coefficient of extinction is negligible.Therefore the correlations necessary for cloud base detection should be calculated for each cloud pixel, p, separately based only on the cloud pixels inside a round moving window centered at each p.This should be done for all pixels that are according to the DEM below p (resulting in ρ below p ) as well as for all pixels at the height of or above p (ρ above p ), resulting in the difference ρ diff p = ρ below p − ρ above p .With ρ above p including p itself, the value of ρ diff p is, according to the above considerations, especially high for pixels that are the lowest that still are located above the local cloud base height.This is because ρ above p is clearly negative here while ρ below p is close to 0 (cf.Fig. 5).Pixels with a local maximum of ρ diff p and clearly negative values of ρ above p can therefore be considered as being most likely cloud base height pixels.
The new algorithm for ground fog detection making use of the theoretical basis developed in this section is presented in detail in Sect.4.1 and is validated in Sect.4.2.As all of the above consideration rely on clouds being cut off by mountains, its application is limited to the Detection Of Ground fog in Mountainous Areas (DOGMA; this acronym will be used hereinafter for the new technique).Since most of the fog occurrence in Taiwan (as well as all the known occurrence of cloud forest) is restricted to the mountain areas that cover the biggest part of the island, the method is suited for ground fog detection in Taiwan despite this limitation.
The optical thickness input and the DEM as well as other inputs that are necessary for the method are described in Sect.3.

Input data and their processing
As shown in Sect.2.2, one of the main inputs for DOGMA is a DEM.The ASTER GDEM 2 (property of METI and NASA), distributed via the USGS global data explorer (United States Geological Survey, 2013) and resampled to a resolution of 250 m, is used.
All other inputs (cf.The cloud optical thickness is based on the MODIS MOD06 2.1 µm optical thickness product derived from radiative transfer calculations (Platnick et al., 2003).As it is only available for daytime scenes, the application of DOGMA is restricted to those (cf.Sect.3.2) Since fog (with the exception of ice fog that can only be found in polar latitudes; Oke, 1978) is completely in the water phase, DOGMA is only applied on water clouds.In addition to the inputs mentioned in Sect.2.2, a cloud mask including information about the cloud phase is thus necessary as an input.It is based on the MODIS MOD35 cloud mask product (cf.Sect.3.3).
To be able to remove pixels that are considerably colder (because they are higher) than the rest of a cloud entity, cloud top temperatures are also needed for the procedure.They are calculated for all water cloud pixels from MODIS thermal infrared (TIR) window channels (cf.Sect.3.4).
To make use of the full 250 m resolution of MODIS, which is only available in the visible bands 1 and 2, different techniques (cf.Sect.3.1-3.4)are applied to sharpen those inputs using the MODIS high-resolution bands.

Pan sharpening of MODIS channels
MODIS channels are not necessary as an input for DOGMA.Some are, however, required to sharpen other MODIS products needed as an input.Therefore several infrared channels are transferred to the 250 m resolution of MODIS channels 1 and 2 using a suitable pan-sharpening technique by Schulz et al. (2012).According to this method a highresolution satellite channel is degraded to match the resolution of the low-resolution channel that is to be sharpened.This is done by averaging each group of 4×4 high-resolution pixels, the area of which corresponds to one and the same low-resolution pixel.Then a potential regression is used to explain the pixel values of the low-resolution channel with the pixel values of the degraded high-resolution channel.The same regression is applied on the original high-resolution channel in order to obtain a high-resolution version of the low-resolution channel.This procedure is, however, not applied on a scene globally as that would lead to a low coefficient of determination of the regression.This would result in a bad quality of the sharpened image.Instead, a channel is sharpened pixel-wise using a different regression for each pixel.Each regression is based on an approximately round moving window with a diameter of 5 pixels centered at the pixel that is to be sharpened.
For DOGMA data pre-processing the pan-sharpening technique has been adapted to MODIS data.Since MODIS has two channels in the resolution of 250 m, both of them were degraded and used to sharpen the low-resolution channels using a multiple potential regression.In order to simulate what the high-resolution image would look like if captured by MODIS in its 1000 m resolution, the sensor's spatial response function is additionally incorporated in the degradation of the 250 m channels.In the averaging process of 4 × 4 high-resolution pixels each high-resolution pixel is weighted by the sensor's spatial response function (taken from Huang The described method is used to sharpen the unitless digital number (DN) values of the MODIS channels 20, 29, 31 and 32 as distributed via the MODIS MOD02/MYD02 product.After sharpening, the DNs of these channels are transferred to radiances using scale and offset values included in the MOD02/MYD02 product.The inverse Planck function is used to calculate black-body temperature (BBT) values from the radiances.
BBTs are also calculated from the 1000 m DNs of the channels 20 and 31.From the DNs of channel 1 the reflectance is calculated using the appropriate scale and offset values from the MOD02/MYD02 product.A sharpening result for channel 31 is exemplarily shown in Fig. 7.An overview of the used channels is given in Table 3.

Cloud optical thickness
Over land surfaces the MOD06 optical thickness is mostly based on channel 1 (Platnick et al., 2003).Therefore a slightly adapted version of the pan-sharpening method used to process the 1000 m MODIS imagery (cf.Sect.3.1) is suited to sharpen it.Instead of a multiple regression a simple regression incorporating only channel 1 (instead of channels 1 and 2) as the independent variable is used.An example result is shown in Fig. 7.

Cloud mask
As it is based on the solar MODIS high-resolution channels 1 and 2 only, the 250 m cloud mask included in the MOD35 cloud mask product does not contain information about the cloud phase.Furthermore it is, for some scenes, heavily flawed in the area of Taiwan.Cities and riverbeds with high reflectance values are often wrongly classified as clouds (cf.Fig. 8).These problems do not occur in the MOD35 1000 m cloud masks.Therefore, the 250 m water cloud mask used as an input for DOGMA is created using the MOD35 1000 m cloud mask as a reference and sharpened as well as unsharpened MODIS imagery as its basis.Before any further processing is applied, thin cirrus are removed from the 1000 m MOD35 cloud mask based on the included cloud classification.Cermak and Bendix (2008) have shown that a threshold applied on the difference BBT TIR − BBT MIR (further referred to as diff 250 and diff 1000 for the 250 m and the 1000 m version) between BBT values in the TIR and medium infrared (MIR) is well suited to distinguish between cloudcontaminated and clear pixels using a threshold.The difference values are clearly negative for cloudy pixels and near to 0 for clear surfaces.For MODIS this difference can be calculated from the channels 31 (TIR) and 20 (MIR).
To identify a suitable threshold for a certain scene several iteratively adapted thresholds are applied on diff 1000 in order to obtain masks in the 1000 m resolution.Those masks are compared to the MOD35 1000 m cloud mask in terms of the percentage of pixels being classified in agreement between the cloud mask derived from diff 1000 and the MODIS reference product.The threshold resulting in the best agreement is chosen.
Since clouds are highly reflective in the visible spectrum, a threshold can also be applied on the reflectance of channel 1 (further referred to as ref 250 and ref 1000 for the 250 m and the 1000 m version) to distinguish between cloud-covered and cloud-free pixels to some degree.This threshold is determined in the same iterative way as the threshold derived from diff 1000 .
These two thresholds can then be applied on diff 250 and ref 250 .The results are two cloud masks in 250 m that still contain several flaws (e.g., cities and river beds being classified as clouds due to their high albedo in the channels 1 and 20).As the flaws often are in different areas the two cloud masks can be combined to a new cloud mask (further referred to as global cloud mask) that consists only of the pixels classified as clouds in both of them.It is of higher quality but still not completely flawless due to the threshold being applied on the full MODIS scene.Therefore thresholds are additionally identified for each 1000 m pixel p separately by using the iterative approach described above for a window of 20 pixels × 20 pixels around p only.The pixels of diff The resulting mask does still not differentiate between water, ice, and mixed phase clouds.Mixed and ice phase clouds are detected using the same threshold approach that is incorporated in the MODIS cloud phase classification.It is mainly based on radiative transfer calculations that have shown that the difference BBT 29 − BBT 31 between the MODIS channels 29 and 31 is high for ice clouds and low for water clouds (Chylek et al., 2006).Also BBTs from channel 31 can be incorporated to detect clouds that are obviously too cold to be in the water phase.The MODIS cloud classification combines these two approaches to identify cloud pixels that are not or not solely in the liquid phase as follows (Platnick et al., 2003): Examples for the enhanced 250 m cloud mask are shown in Fig. 8.

Cloud top temperature
Cloud top temperatures are calculated for all water cloud pixels from the sharpened BBTs of the MODIS TIR window channels 31 and 32 using a split-window algorithm proposed by Jiménez-Munõz and Sobrino (2008).While the method was not explicitly developed for cloud surfaces, there is no reason why it could not be applied here.It allows us to correct for atmospheric absorption as well as emissivity effects.The latter, however, are ignored in our cloud top temperature retrieval as an emissivity of 1 for the cloud surface is assumed.This approximation for water clouds is possible for TIR wavelengths, if the cloud is supposed to have a thickness of at least several tens of meters (cf., e.g., Yamamoto et al., 1970;Hunt, 1973).The total atmospheric water vapor content, which is needed as an additional input for the cloud top temperature retrieval, is taken from the MOD05/MYD05 total precipitable water product.It is resampled to the 250 m resolution without sharpening.

DOGMA -detailed description
DOGMA is fed with the inputs described in Sect. 3 (cf.Fig. 9).Pixels classified as ice cloud or mixed phase cloud in the 250 m cloud mask are removed from the analysis and marked as unclassifiable as they might block the view to lower fog layers.The algorithm runs through each 250 m water cloud pixel p of a MODIS scene and calculates ρ below p , ρ above p , and ρ diff p (cf. Fig. 10a) from a round window with a diameter of 40 pixels around p as described in Sect.2.2.The scene is then scanned for local maxima of ρ diff p in order to detect pixels at the cloud base height (further referred to as CBH pixels).This is done by comparing each water cloud pixel p to all other water cloud pixels in a round window with a diameter of 20 pixels centered at p. Pixels which are higher than the lowest direct neighbor of p or lower than the highest direct neighbor of p are excluded from this comparison as they would (if p were actually a CBH pixel) most probably be CBH pixels of the same cloud base as p. p is marked as a low-certainty CBH pixel (cf.Fig. 10b) if a. p has a higher ρ diff p than all the pixels it is compared to b. ρ diff p is above 0 c. ρ above p is clearly negative (an empirically derived threshold of −0.3 is used) and d. p has, according to the DEM, a slope elevation of at least 7.2 % (this is checked as the considerations in Sect.2.2 are all based on the assumption of mountainous terrain).
As ρ diff p and ρ above p are calculated from relatively small windows, small-scale variations of the cloud bottom height can be captured well.On such a local level random smallscale gradients of the cloud thickness or the extinction coefficient that spatially coincide with an increase or decrease in terrain height could be mistaken for correlations that are caused by a cloud being cut off by the terrain.For bigger windows, resulting in a bigger sample size of the correlations, such small-scale gradients have a much lower impact.Therefore ρ above p is calculated again (further referred to as ρ above p, 120 ) for each low-certainty CBH pixel p.This is done based on the optical thickness and height of all pixels in a round window with a diameter of 120 pixels centered around p. For a window of this size the assumption of a negligible variability of cloud top height, cloud bottom height and coefficient of extinction inside the window may not be fulfilled.If p is actually a CBH pixel, the 120 pixels window could, for example, contain water cloud pixels with a higher elevation than p that are free of ground fog due to variations of the cloud base height.Therefore the correlation may not reach very low values even for CBH pixel.Thus, a relatively high threshold of 0 is applied on ρ above p, 120 .If ρ above p, 120 is below that value, p can be regarded as a CBH pixel with medium certainty.
As real CBH pixels should be a part of a cloud base that is formed by several CBH pixels, CBH pixels that are surrounded by other CBH pixels should have an increased probability of being real CBH pixels.Therefore for each mediumcertainty CBH pixel p, the algorithm checks whether at least 10 other medium-certainty CBH pixel can be found in a round window with a diameter of 40 pixels centered at p.If that is the case, p is marked as a high-certainty CBH pixel.
The filtering that is necessary to obtain high-certainty CBH pixels is often too strict, resulting in too many pixels being filtered out.In order to increase the number of detected CBH pixels (which is necessary to capture small-scale variations in the cloud bottom height as precise as possible),  DOGMA goes back to the original low-certainty CBH pixels and filters them again based on their vertical distance to an assumed cloud base surface that can be defined from the high-certainty CBH pixels.This cloud base surface is modeled for each water cloud entity separately by interpolating the heights of all high-certainty CBH pixels belonging to that cloud entity using inverse distance weighting (IDW; Shepard, 1968).All low-certainty CBH pixels that are within a vertical distance of less than 400 m to this surface (further referred to as final CBH pixels) are used for the discrimination of ground fog.
For the final ground fog discrimination, the heights of all final CBH pixels of each water cloud entity as well as their temperature taken from the cloud top temperature image are interpolated using IDW interpolation (cf.Fig. 10c).
Each pixel p in which the interpolated height of the cloud base is below or equal to the height taken from the DEM is a ground fog pixel if the IDW interpolated temperature of p is not more than 3 K higher than the temperature taken from the cloud top temperature image.If the measured temperature was much colder, p would be significantly higher than the pixels used for the IDW interpolation.Therefore, it would either belong to another cloud level or the assumption of a cloud top surface that can locally be approximated as plane would be wrong for p.
So far the existence of a cloud base surface that has a onedimensional intersection with the terrain has been assumed.If a fog cloud does, however, touch the terrain with its entire base (complete valley fill), such an intersection does not exist.Therefore the cloud would not be identified as fog by the tests described so far.In order to detect those clouds nevertheless, each cloud entity from the water cloud mask is further examined when no fog has been detected in it.For each pixel p of the entity, ρ is calculated from the pixels in a round window with a diameter of 40 pixels that is centered at p.If the median of all ρ is below −0.3, all pixels of the entity are considered ground fog pixels.The final ground fog mask is shown in Fig. 10d.

Methodology of validation
For many areas, METeorological Aerodrome Reports (METARs) are ideally suited to validate ground fog detection schemes (cf., e.g., Cermak and Bendix, 2011;Schulz et al., 2012).These reports are based on weather observations from airports and include cloud base heights as well as information about the visibility at ground level.In Taiwan, however, all airports are located in the plains in the outer parts of the island and not in the mountainous area where DOGMA can be applied (dark red area Fig. 10b).
Considering this lack of data, several cameras of the type PlotWatcher Pro (Day 6 Outdoors, LLC, USA) were installed in the mountainous parts of the island.The cameras are set up to take a photo each minute and save it to a SD card.Additionally, the data of a CL31 ceilometer (Vaisala, Finland) installed in Xitou have been used (cf.Fig. 11; Table 4).Table 4 gives an overview about the time spans during which the different instruments have captured data.For those time spans DOGMA ground fog products were generated from the data of all Terra and Aqua daytime overflights over Taiwan.The DOGMA ground fog product (cf.Sect.4.2.1) as well as the IDW interpolated cloud base surface (cf.Sect.4.2.2) have been compared to the camera and ceilometer data in order to assess their quality.

Validation of the ground fog product
The cameras 1 to 6 are located at positions that are often cloud immersed.All of their images taken at the time of MODIS overflights were manually classified into the categories "fog free" or "fog immersed".From the ceilometer data captured at overflight times, information about the cloud immersion of the instrument was also extracted.Each camera and ceilometer reference observation was then compared to the DOGMA ground fog product pixel that is corre- spondent in space and time.Ice cloud and mixed phase pixels have not been included in this comparison as DOGMA marks them as unclassifiable.The results have been summarized in a contingency table.From this table the following statistical measures have been calculated (cf.Mason (2003); Matthew (1975)

Validation of the DOGMA cloud base height
As the quality of the DOGMA ground fog product is a direct result of the cloud base surface height (cf.Fig. 10c), the latter has been validated separately.This was done using the ceilometer as well as camera 7. The ceilometer is located in the valley of Xitou where a high fog frequency can be observed.The camera is located near the lower end of a valley on Chi-Lang Mountain facing up the valley.While the camera location is usually fog free, the upper parts of the valley are often fog immersed.The intersection between the cloud base and the terrain can be observed by the camera in these cases.

Cloud base validation using ceilometer data
The ceilometer validation has been conducted for all MODIS overflights in which the DOGMA cloud base product provides information for the ceilometer location and the ceilometer provides cloud base information (this means that it is not cloud immersed and the pixel it is located in is not cloud free).For these overflights, the DOGMA cloud base height has been extracted for the pixel of the ceilometer location.The mean deviation has been calculated between ceilometer reference data and the extracted cloud base heights.All scenes in which the ceilometer obtained cloud base height is above the height of the highest pixel in the valley of Xitou that is cloud covered (according to the 250 m water cloud mask) were excluded from this calculation.This was necessary as DOGMA extracts the height of the cloud base from the DEM and is therefore -as a matter of principle -not able to detect any cloud base that does not touch the terrain.For ground fog detection information about higher cloud bases is not of interest.Scenes in which the DOGMA cloud base is below the ceilometer location also needed to be excluded from the calculation as ceilometer data (which are only available when the cloud base is above the ceilometer) would only be available if the DOGMA cloud base height were wrong.This would result in a biased validation.

Cloud base validation using camera data
While cloud bottom heights recorded by the ceilometer could be directly compared to the height of the DOGMA cloud base, the camera footage needed to be processed manually in order to obtain reference cloud bottom heights.All images taken at MODIS overflights were assessed to determine whether the captured slopes of Chi-Lan Mountain are cloud immersed, as shown in Fig. 12. Scenes for which that is not the case and scenes for which the weather conditions (e.g., ground fog at the camera location) did not allow an unequivocal identification of the cloud base were excluded from further analysis.For the remaining scenes the intersection of the cloud base height and the terrain was marked manually in the image (red lines in Fig. 12).For those image pixels marked as intersection pixels, the altitude was extracted from the ASTER GDEM 2 reprojected to the view of the camera (cf.Schulz et al. (2014) for details about the reprojection).IDW interpolation has been used to model a cloud base surface from the heights and positions of those pixels in a resolution of 250 m for a domain defined by their bounding box.For each MODIS scene for which the cloud base obtained from the camera data has been calculated and the DOGMA cloud base product also contains information inside the bounding box domain, the median of the deviations in height between all pixels of the camera obtained cloud base surface and the corresponding DOGMA cloud base surface pixels was then calculated.The result is a single value for each scene that describes the deviation between the camera obtained cloud base and the DOGMA cloud base for the whole view shed of the camera.From those deviations the mean deviation has been calculated.

Validation results and discussion
The results of the validation of the DOGMA ground fog product are shown in Tables 5 and 6.Also the results of a validation carried out by Cermak and Bendix (2011) for their own method using METAR data from a European domain are included.As the two validations are based on different data sets from different areas, they should not be directly compared to each other.The validation by Cermak and Bendix should only be seen as a reference for the current quality of ground fog detection from spaceborne sensors.
As shown in Table 6 both methods have a relatively similar overall fog detection quality (MCC).DOGMA tends to underestimate the fog frequency while the method by Cermak and Bendix generally overestimates it (bias).In detailthat means not only that a lower ratio of fog pixels is classified correctly as fog by DOGMA than by the method of Cermak and Bendix (POD) but also that a lower ratio of all pixels that are classified as fog contaminated is wrongly classified (FAR).A higher ratio of fog-free pixels is, however, wrongly classified as fog by DOGMA (POFD).The main issue with both methods is a relatively high FAR in combination with a relatively low POD.
The shortcomings of DOGMA might mostly be caused by clouds that cannot locally be approximated as planeparallel and horizontally homogeneous as it is assumed by the method.Also, clouds with a very high optical thickness that often cover Taiwan can cause problems.The footage of the cameras shows that optically thick clouds do often touch the ground, even if it is not possible for a human observer to distinguish ground fog in a satellite image.Also, DOGMA has problems detecting fog in these situation.An optical thickness of 40 -a value that is surpassed in 13.31 % of the pixels that have been used for the validation -corresponds to a transmittance of 4.24 × 10 18 .With such a low transmittance almost no light -and therefore hardly any information -from the cloud base reaches the sensor.In order to assess the impact of those very thick clouds, a validation excluding all reference observations during which the camera or ceilometer has been located below a cloud with an optical thickness of more than 40 has been carried out.As shown in Table 6, this increases the overall detection quality (MCC), mainly due to a better POD.
The mean deviation of the DOGMA cloud base height is given in Table 7. DOGMA does not directly calculate a cloud base height.Instead CBH pixels are detected in a two-dimensional image (cf.Sect.4.1).As the validated DOGMA cloud base height is interpolated from the height of the DEM extracted from these pixels, its precision is (if inaccuracies of the interpolation are ignored) a result of the correctness of the detection of CBH pixels in combination with the topography and its mapping in the 250 m resolution of the DEM.The mean steepness of the slopes of both valleys incorporated in the cloud base validation is approximately 50-60 %.As a pixel with sides of 250 m in length has a diagonal of about 354 m, height differences of 177 m (50 % of 354 m) to 212.4 m (60 % of 354 m) inside a single pixel are possible and height differences between 125 m (50 % of 250 m) an 150 m (60 % of 250 m) are unavoidable.This means that even CBH pixels that have been detected in the correct position may result in relatively imprecise cloud base heights, although the delimitation of the fog immersed area would be perfect.Conversely, this implies that the ∼ 200 m error of the DOGMA cloud base is the result of a relatively small mean error of less than two pixels in the horizontal positions of the CBH pixels.Therefore the DOGMA cloud base height product is suited for ground fog delimitation but should not be used for other purposes.

Conclusion and outlook
Common fog detection schemes that have been developed for radiation fog are not applicable in Taiwan, as they rely on assumptions (cf.Sect.2.1) that are not met by most fog occurrences in Taiwan.Therefore the presented method has been developed.DOGMA does not calculate the cloud base height from the difference between the cloud top height and the cloud thickness.Instead, pixels at the cloud base are directly detected using a statistical approach that is based on a negative Spearman's rank correlation coefficient between the optical thickness and the terrain height of fog immersed pixels extracted from the DEM.As it relies on clouds that can be at least locally approximated as plane-parallel and horizontally homogeneous, and the MOD 06 optical thick- ness product is based on radiative transfer calculations using a plane-parallel cloud model (Platnick et al., 2003), DOGMA does rely on some assumptions.The necessary degree of plane-parallelism and horizontal homogeneity is, however, low enough so that the method is applicable to seas of clouds (cf.Fig. 2), which are typical in Taiwan.As the comparison to the method of Cermak and Bendix (cf. Sect. 5) has shown, the overall quality of DOGMA's fog detection is (despite some problems with fog clouds with a very high optical thickness) comparable to that of a modern fog detection scheme developed and validated for the temperate zones.It therefore seems to be applicable to the creation of ground fog frequency maps that can be used for the country-wide mapping of Taiwan's cloud forests.The relationship between fog frequencies and the occurrence of mountain cloud forest in Taiwan will be the subject of future work.In addition, it will be investigated whether a daytime-only approach is suited for the delimitation of mountain cloud forest.
DOGMA is restricted to mountainous areas.If perfectly plane-parallel and horizontally homogeneous fog clouds are assumed, it should work even for slightest slopes.As those perfect conditions would not actually be met, it will be the subject of future work to find out for which areas DOGMA is suited.It would, for example, be conceivable that the method works for very plain radiation fog in the valleys of low mountain ranges.
The vertical precision of the DOGMA cloud base is highly dependent on the horizontal accuracy of the fog detection (cf.Sect.6).Also, fine details in the margin of fog entities that are restricted in their extent by the terrain can be much better captured in the highest MODIS resolution of 250 m than in the 1000 m resolution (cf.Fig. 7a).For those reasons all inputs necessary for DOGMA needed to be transferred to the 250 m resolution.An inclusion of a 250 m version of the optical cloud thickness product in future MODIS data collections could make such a step superfluous and would be useful for future ground fog detection schemes.
As the only necessary inputs for DOGMA are a DEM, a cloud mask, an optical thickness product, and cloud top temperatures, its implementations for the imagery of other satellites than MODIS should be possible without problems.MODIS has been chosen for the first implementation of DOGMA due to its high spatial resolution and the long time span for which its data are available.The drawback of a polar-orbiting satellite is, however, its bad temporal coverage.Therefore the imagery of the new Himawari 8 satellite with a resolution of up to 500 m pixel −1 and a sampling rate of 10 min would be well suited for future implementations of DOGMA in order to obtain information about ground fog in Taiwan.

Figure 1 .
Figure 1.Ground fog detection under the assumption of a planeparallel cloud geometry.

Figure 2 .
Figure 2. Surface temperatures calculated using a split-window approach by Jiménez-Munõz and Sobrino (2008) for the MODIS overflight on 5 January 2014 at 10:35 UTC + 8.The emissivity for the land surface has been taken from the MODIS MOD 11 product.For clouds an emissivity of 1 has been assumed (cf.Sect.3.4).

Figure 3 .
Figure 3. Theoretical geometrical thickness/optical thickness under the assumption of a perfectly plane-parallel cloud restricted by the terrain in its extent.The assumed cloud reaches from a base height of 900 m a.s.l. up to a top height of 1400 m a.s.l.

Figure 4 .
Figure 4. Relationship between the optical thickness of a cloud partially touching the ground and the terrain height in Xitou, Taiwan, on 31 November 2014.The lower plot is based on all cloud pixels inside the plot domain.For pixels above the ceilometer-derived cloud base ground contact of the cloud can be assumed under the assumption of plane-parallelism.The cloud base height was detected by a ground-based ceilometer (cf.Sect.4.2 for information about the ceilometer and the location).The ASTER GDEM 2 resampled to a resolution of 250 m (cf.Sect.3) has been used for the terrain height.The optical thickness was taken from a sharpened 250 m version of the MODIS MOD 06 product (cf.Sect.3.2).

Figure 5 .
Figure 5.The calculation of ρ diff p for pixels p in different locations.(a) Inside a cloud, ρ above p and ρ below p are calculated from pixels with ground fog.Therefore both values are clearly negative.This results in a small ρ diff p .(b) At the cloud base height, ρ above p is calculated from pixels with ground fog.It is therefore clearly negative.ρ below p is calculated from pixels without ground fog.Its value is therefore close to 0. This results in a high ρ diff p .(c) Below a cloud, ρ above p and ρ below p are calculated from pixels without ground fog.Therefore both values are near to 0. This results in a small ρ diff p .
Fig. 6) are based on MODIS Collection 051 Level 1B and Level 2 products.MODIS data have been chosen instead of the data of geostationary satellites (e.g., the Japanese Himawari series covering the area of Taiwan) because of the combination the long time span for which the data are available (MODIS: since 1999 (Terra)/2002 (Aqua); Himawari 7: since 2006; Himawari 8: since 2015) and their relatively high spatial resolution (MODIS: up to 250 m; Himawari 7: up to 1000 m; Himawari 8: up to 500 m).The latter is necessary due to the complex topography of Taiwan's mountains.

Figure 6 .
Figure 6.Flowchart for the processing of the inputs required to run DOGMA.
250 and ref 250 covered by p are then classified as cloud contaminated or cloud free using these thresholds.The results from the local approach and the global cloud mask are combined to a final 250 m cloud mask.A pixel is considered as cloudy if clouds are present according to the global cloud mask and a. the pixel is cloud contaminated according to the local thresholds applied on diff 250 and ref 250 or b. the pixel is unambiguously cloudy according to the value of diff 250 (diff 250 is below the diff 250 value of at least 50 % of the pixels in the 20 pixels × 20 pixels window that are considered as clouds according the MODIS 1000 m cloud mask).

Figure 10 .
Figure 10.Steps in the calculation of the DOGMA fog mask for the same MODIS overflight as shown in Fig. 2 (5 January 2014, 10:35 UTC + 8).(a) ρ diff p calculated for each pixel of the 250 m cloud mask.(b) Low-certainty CBH pixels.The dark red area is not steep enough to apply DOGMA.(c) IDW interpolated cloud base height.(d) Final DOGMA ground fog mask.

Figure 11 .
Figure 11.Locations of the instruments used in the validation study.

Figure 12 .
Figure 12.Manual identification of cloud base heights on Chi-Lan Mountain using PlotWatcher Pro imagery (scene from 31 August 2013).

Table 1 .
Different cloud height retrieval methods.

Table 2 .
Different cloud thickness retrieval methods.

Table 3 .
MODIS channels that are used in the creation of input data for DOGMA.
et al., 2002)of the low-resolution pixel that is to be sharpened.

Table 4 .
Overview of the instruments used for the validation study.

Table 5 .
Confusion matrices for the validation of the DOGMA ground fog product.

Table 6 .
Validation results for the DOGMA ground fog product.

Table 7 .
Validation results for the DOGMA cloud base height.