Topoclimate Mapping Using Landsat ETM+ Thermal Data: Wolin Island, Poland

: Variations in climatic pattern due to boundary layer processes at the topoclimatic scale are critical for ecosystems and human activity, including agriculture, fruit harvesting, and animal husbandry. Here, a new method for topoclimate mapping based on land surface temperature (LST) computed from the brightness temperature of Landsat ETM+ thermal bands (band6) is presented. The study was conducted in a coastal lowland area with glacial landforms (Wolin Island). The method presented is universal for various areas, and is based on freely available remote sensing data. The topoclimatic typology obtained was compared to the classical one based on meteorological data. It was proven to show a good sensitivity to changes in topoclimatic conditions (demonstrated by changes in LST distribution) even in flat, agricultural areas with only small variations in topog ‐ raphy. The technique will hopefully prove to be a convenient and relatively fast tool that can im ‐ prove the topoclimatic classification of other regions. It could be applied by local authorities and farmer associations for optimizing agricultural production. Contributions: Conceptualization, M.C.; methodology, J.C.; tion, analysis, M.C.; investigation, resources, curation, M.C.; writing— preparation, M.C.; writing—review editing, J.C.; visualization, M.C.; supervision, M.C.; administration, M.C.; funding


Introduction
Local climates are controlled by global, regional, and local factors. For example, planetary waves operate in large to macro scales (>10 3 km) and high-pressure systems or trajectories of cyclones contribute to mesoscale processes (10 1 -10 3 km). The most common variation in climatic patterns, however, is at topoclimatic scales between 10 1 and 10 −3 km due to boundary layer processes. The boundary layer is the part of the lower troposphere that interacts directly with the earth's surface through turbulent transport processes. Topoclimate, meaning local climate, is controlled by such factors as relief, vegetation, hydrographic conditions, and soil. On a larger scale, these local factors such as vegetation type cease to be climatic factors and become phenomena dependent on the climate.
How climate factors are constructed has a significant effect on species distribution models and climate-change predictions. Topoclimates may cause significant variations (p value < 0.002) in the distribution of ferns or grasses that cannot be explained by macroclimatic variables [1]. Models where factors in topoclimatic scale such as cold air drainage, topography, or habitat are not considered in the climate mapping methodology may result in misleading conclusions [1]. Therefore, [2] suggested that future studies on species distribution models should use the temperature dataset in topoclimatic scale that best reflects the ecology of the species, rather than automatically using low-resolution data from WorldClim (0.9-18.6 km by pixel) or CHELSA (~1 km/px).
TRS allows monitoring of the thermal properties of the earth's surface [20]. Thermal satellite sensors range from low resolution (500-1000 m by pixel), e.g., Sentinel 3 Sea and Land Surface Temperature Radiometer (SLSTR) of 500 m/px, Moderate Resolution Imaging Spectroradiometer (MODIS) and NOAA Advanced Very-High-Resolution Radiometer (AVHRR) of 1000 m/px), to medium resolution e.g., Landsat Thematic Mapper (TM)-120 m/px, Landsat 7 with Enhanced Thematic Mapper Plus (ETM+)-60 m/px, Thermal Infrared Sensor (TIRS) of Landsat 8-100 m/px, and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)-90 m/px [21][22][23][24][25][26][27]. As the accurate LST estimation requires satellite data of a good spatial resolution satellite data, Landsat 7 thermal band 6 (10.4-12.5 μm) with the best spatial resolution of 60 m in thermal bands was used. Landsat 7 images consist of eight spectral bands with a spatial resolution of 30 m for Bands 1-5 and 7. The resolution for Band 8 (panchromatic) is 15 m. Whereas most bands can collect one of two gain settings (high or low) for increased radiometric sensitivity and dynamic range, thermal data from band 6 is acquired in two bands from one detector in both high (Band 6H) and low (Band 6L) gain for all scenes, which results in a spatial resolution of 60 m/px. We have used CLC2006 to match with Landsat ETM+ data.
This research was performed in a coastal area, Wolin Island, Poland, which features a temperate climate. The results were verified through comparison with Paszynski's classification [28], which is broadly used in Poland. A digital version was elaborated by us for the purposes of this study.

Study Area
Wolin, the largest island (265 km 2 ) of Poland, is located between the Szczecin Lagoon, the mouth of Odra River, and Pomeranian Bay ( Figure 1). The geographical uniqueness of this region lies in its landscape diversity, insular character, and diverse flora and fauna [37]. The LULC form map (Figure 2) of Wolin Island shows that coniferous forests occupy the largest area of the island (25%), mostly in its northeastern and western parts. Nonirrigated arable lands (17%) and pastures (14%) occur in the eastern part of the island. Mixed forests (12%) and broadleaf forests (8%) prevail in the northern part of the island and the Świna Gate region (Figure 1) in the west. Inland marshes (7%) occur mainly in the Rów Peninsula, Świna Gate, and Dziwna Spit. Complex cultivation (3%) and land principally occupied by agriculture (including significant areas of natural vegetation) (3%) are dispersed in various locations across the island. Discontinuous urban fabric represents 2% of the island's area, whereas other LULC forms are <1%. The steepest slopes of the island (up to 54°) are related to moraine hills and coastal cliffs ( Figure 3). Gentle slopes and flat areas predominate, mostly in the southeastern and western part of the island, forming a 1.7° average slope for Wolin Island. Overall, flat areas (slopes of <2°) occupy 74% of the island's total area, gentle slopes (2-5°) occupy 16%, and steep slopes (>5°) occupy only 10%. Slopes with eastern (32%) and western (29%) exposures prevail over the slopes with southern (21%) and northern (18%) exposures. According to the Köppen-Geiger climate classification [38][39][40], Wolin Island's climate is categorized as Cfb, representing a warm temperature and year-round precipitation (no dry season, warm summer).

Data and Methods
One of the main goals of this research was to present a technique that will become a convenient, relatively easy, and fast tool to improve topoclimatic classification that can be applied by a wide group of users. Therefore, we developed the method to use free datasets and free software, if available (including CORINE Land Cover and Landsat images, and free software: SAGA GIS, LST program, and Map Comparison Kit).

Selection of Landsat Images
Landsat Enhanced Thematic Mapper Plus (ETM+) images were used for the LST map calculations ( Figure 4). The data were acquired from the United States Geological Survey (USGS) Resource Center (http://glovis.usgs.gov, accessed on 26 June 2021). ETM+ is the main instrument aboard Landsat 7 and it records thermal infrared images with resolutions of 15 (panchromatic), 30 (multispectral), or 60 m per pixel. Some images are not suitable for the analyses, including those with (1) too high a percentage of cloud cover or (2) that are outdated with no reference to the relatively current map of LULC used in the research (in terms of the distribution of settlements or road networks on the island, for example). Moreover, only images taken before May 2003 were used because of the failure of the Scan Line Corrector (SLC) instrument on 31 May 2003 (http://landsathandbook.gsfc.nasa.gov/payload/prog_sect3_3.html, accessed on 26 June 2021). Images taken at different times of the year (from March to September, Figure 5) were selected to increase the universality of the typology created. Altogether, seven Landsat ETM+ images taken between 9:50 and 9:55 GMT were selected (Table 1).   [19] with the algorithm of Qin et al. [42], making use of the Landsat ETM+ thermal bands. All the images are presented using the same temperature scale (−3-55 °C) so that the spatial variability of the LST values are easier to compare.

Retrieval Algorithm
LST from Landsat thermal band 6 data is commonly retrieved through the algorithm proposed by Qin et al. [42]. The algorithm is based on the linearization of Planck's radiance function and is expressed as follows: where Ts is the land surface temperature in K, Tsensor is the brightness temperature in K computed from Landsat ETM+ band 6, Ta is the effective mean atmospheric temperature (K), and a6 and b6 are constants with values of −67.355351 for a6 and +0.458606 for b6 when the LST is between 273.5 and 343.5 K. C6 and D6 are calculated using the following equations: where ε is the ground surface emissivity and τ6 is the atmospheric transmittance ( Table 2). Following Qin et al. (2001), τ6 is estimated from the atmospheric water content (Table 3).  Table 3. Ratio of the distribution of water vapor content to the total content of water vapor as a function of altitude (Rw(z)) in different four standard atmospheric profiles [42]. Ta is calculated from one of linear Equations (4)-(7) referring to the four standard atmospheres:

Altitude
Ta = 16.0110 + 0.92621T0 (for mid-latitude summer), where T0 is the near-surface air temperature. Temperatures measured at meteorological stations can be used. In this study, Equations (6) and (7) were used for mid-latitudes. The total atmospheric water vapor content (w, g cm −2 ) is calculated as follows: , where w(z) is the water vapor content at altitude z. Rw(z) is the ratio of water vapor content at a given altitude to the average for a standard atmospheric column and is given in Table  4. Water vapor at altitude z (absolute humidity) was calculated from the formula 219.5 (9) where g is the absolute air humidity (g m −3 ), e is the actual water vapor pressure (hPa), T is obtained from meteorological data, and e is calculated by converting the formula for relative humidity [14]: 100%. Therefore, where f is the relative humidity (%) and emax is the maximum water vapor pressure (%), which is read from the psychometric boards. Ground surface emissivity (ε) was calculated using the LST software [19], which has a built-in emissivity calculation module. Since plants and trees control their own temperature through evapotranspiration, the normalized difference vegetation index (NDVI) was introduced to the applied module. The LST calculation based on the normalized difference vegetation index (NDVI) from Landsat ETM+ is expressed as follows: The relation between the emissivity and the NDVI (for NDVI values from 0.157 to 0.727) can be expressed by the following equation: 1.0094 0.047 ln NDVI .
For other NDVI values, fixed emissivity values were assigned. For most of the island area, the NDVI values were in the range 0.157-0.727. Only areas of water exceeded this range, and the value of 0.985 was assigned to water in accordance with the work of Snyder et al. [43].
The LST maps are presented in Figure 5. Forested areas are characterized by lower LST values in the mornings, since vegetation isolates the ground from incoming radiation. Water is characterized by the lowest LST values (Figures 1 and 5), except during August and September, when the water temperature is close to the temperature of forested areas due to warm water in the summer. Urban areas are characterized by the highest LST, especially in July. Agricultural areas also have higher LST values, especially in August and September; the bare soil heats quickly after crops are harvested and grasslands are grazed or mowed.
All images were taken during sunny and stable weather conditions lasting for at least several days. According to the Grosswetterlangen classification of atmospheric circulation patterns [41], the meridional circulation group occurred five times (in five Landsat ETM+ images) and the mixed circulation group occurred two times. However, the seven Landsat images represent six different circulation types (Table 1) within these two groups. This is important for the universality of the approach: despite the limited number of images used in the analysis, they represent different seasons of the year, different weather conditions, and are statistically different (Table 4)-a linear correlation between each pair of images (datasets) is weak. This makes our conclusions more generalizable.

Unsupervised Classification of LST Images
The classification of thermal images was performed in order to find areas with similar LST values. This enabled the indirect delimitation of the topoclimates. The unsupervised classification was performed on standardized data and relied on mathematical differences between LST values-the concept was to find LST distribution patterns based on computer learning. It was not possible to arbitrarily establish the thresholds for the topoclimatic classes, as the differences between them were difficult to establish on the basis of visual analysis only. The classification obtained was compared to previously published topoclimate distribution based on the other approach [28,44]-see Section 3.5.
The classification was performed based on seven separate thermal images (not averaged). The standardization function of the STATISTICA software was used to standardize the data. Unsupervised classification was initially carried out with the use of the isoclust, isodata, and cluster [45] procedures.
After the unsupervised classification process, the significance of the LST differences within designated classes was tested with a variance analysis, post hoc tests, and discriminant analysis, as described below. First, a one-way analysis of variance (ANOVA) was used. A low value of testing probability (p-value) indicates the significant dependency of the dependent variables (LST images) tested on the independent variable (classification map). The analysis of classification maps showed that the p-value was close to 0 (p = 1 × 10 −6 ), which means that classes differ significantly in terms of thermal properties (Table 5). The second tool, post hoc tests, informs on which LST classes differ significantly. The post hoc Scheffe test in the STATISTICA program [46] was used. With few exceptions, this test showed that the classes differed significantly. For example, on the ISOCLUST8 map with seven classes, only 1 combination from 343 (i.e., 7 images out of 49 combinations) was not significantly different (Table 6). The third tool, discriminant functions, was applied to determine whether the variables discriminated between two or more naturally occurring groups. The Wilks' lambda discriminant function for all seven images was analyzed. Wilks' lambda is a standard statistic used to denote the statistical significance of a model's discriminatory power. Wilks' lambda assumes values from 0 (perfect discrimination) to 1 (no discrimination). The isodata classification method ( Figure 6) gave the most coherent and unambiguous results, with a low Wilks' lambda value and a high percent of clearly classified pixels compared to other methods (Table 7). Therefore, it was chosen for use in further analyses. Figure 6. The Wolin Island area LST classification map (ISODATA11) using isodata and the k-means procedures. Isodata and k-means methods are described in [47] and [48].
The LULC form definitions were based on the CORINE Land Cover 2006 (CLC2006) database (http://www.eea.europa.eu/publications/COR0-part1, accessed on 26 June 2021; Figure 2). The description of classes is generalized, and only LULC forms that occupy > 10% of the class were included in the description. Therefore, the percentages of various LULC forms do not necessarily add to 100%.
Thermal type II, which is 60% agricultural terrain, was named "built-up areas": small villages were classified in the CLC2006 database as agricultural areas/pastures even when they lie next to agricultural areas/pastures; this is due to the limited accuracy of the CLC2006 database. The authors separated villages based on the visual interpretation of the aerial and motor-glider low-altitude visible images [49], and topographic maps that show the spatial distribution of urban fabric on Wolin Island.
The lithology/moisture map (Figure 7) is the combination of a lithology map and the moisture map, which consists of 11 classes (Table 8). The lithology was adapted in a generalized form by reducing the number of classes from the detailed geological map of Poland at a scale of 1:50,000, prepared as single map sheets [50][51][52][53][54][55]. Only the main types of land cover that differed significantly in terms of thermal properties [56][57][58], such as specific heat, heat capacity, thermal conductivity, and thermal diffusivity, were kept for further analyses (Figure 7).  Since a moisture map of Wolin Island is not available, a 1 m hydroisobath was digitized from the hydrographic map of Poland at a scale of 1:50,000 [59]. Moisture information is critical, as the water content determines the thermal properties of the ground due to the specific thermal properties of the water. It was assumed that the areas with shallow groundwater tables are generally wet, while the areas with a deep groundwater table are generally dry.

Digital Elaboration of Paszyński's Classification
The method of determining topoclimatic maps developed by [44] is based on the characteristics of energy exchange between the atmosphere and the ground. This is a fundamental process that shapes local climates. The Paszyński classes were originally separated on the ground as a result of direct heat flux measurements, and then described based on relief (slope and aspect) and land cover characteristics (Figure 8 and Table 10). Energy exchange at the interface between the atmosphere and the ground can be presented quantitatively in the form of the energy balance equation, also known as the heat balance equation. The simplified equation is Q * ± H ± E ± G = 0, (14) where Q * is the radiation balance, H is the turbulent perceived heat flux to or from the atmosphere, E is the latent heat flux associated with the processes of evaporation or condensation, and G is the heat flux in the ground (mainly conducted). Component signs (+ or −) depend on the direction of the flux from the atmosphere to the active surface or vice versa. The final topoclimate classification, taking into account day and night, includes 16 major topoclimatic types [44] corresponding to the relative values of the major components of heat balance. This classification refers essentially to the growing season and stable, sunny weather. Each of the 16 topoclimate types were also described by Paszyński in terms of relief and land cover/land use characteristics. Paszynski's division was elaborated and digitalized in this report to verify the LSTbased topoclimate categorization method. Relief characteristics (slope and aspect) describing topoclimates types were calculated from DTED 2 and then reclassified in accordance with Paszyński's description. Concave forms were defined based on the "r.geomorphon" extension [62] implemented in the GRASS GIS software [63].
LULC classes were designated based on the CLC2006 database, which distinguishes classes within Wolin Island. However, we simplified Paszyński's topoclimatic classification [28,44] by eliminating or merging classes with minor relevance or ambiguous definition-e.g., class 1.4 (areas varied in terms of topography) was not included in the digital elaboration. In total, 13 topoclimate classes were thus distinguished for Wolin Island (see Figure 8, Table 10). Water 1 Figure 8. Terrain topoclimatic classification according to [28,44] and based on the digital elaboration of the relief and LULC forms. 1: Water; 2: south exposition and slope >5; 3: all areas exclusive north and south exposition and slope > 5°; 4: north exposition and slope > 5°; 5: flat forms, slopes < 5°; 6: concave forms; 7: south exposition, slope > 5°; 8: flat forms; 9: all areas with exclusive north and south exposition, slope > 5°; 10: north exposition, slope > 5°; 11: urban and industrial areas with convex forms; 12: urban and industrial areas with plains; 13: urban and industrial areas with valleys.

Map Comparison Tools
Cross tabulation was used to compare the two derived typologies, LST map, and the modified Paszyński's map described in Section 3.5. Cross tabulation allows for the evaluation of similarities between two maps. The kappa index of agreement (KIA) [64] was calculated. It indicates the degree of agreement between two maps. KIA itself is usually not sufficient for a classification scheme that attempts to accurately specify both quantity and location. The VALIDATE module (IDRISI), which was also used, offers one comprehensive statistical analysis that describes how well a pair of maps agree in terms of (1) the quantity of cells in each category and (2) the location of cells in each category. To make this comparison, it offers three options: kappa for no information (Kno), kappa for gridcell level location (Klocation), and kappa for stratum-level location (KlocationStrata). Moreover, VALIDATE divides agreement and disagreement into the more sophisticated components as agreement due to chance, agreement due to quantity, agreement due to location at the stratified level, agreement due to location at the grid cell level, disagreement due to location at the grid cell level, disagreement due to location at the stratified level, and disagreement due to quantity [31].
For a better understanding of agreement and disagreement of quantity and location, assume that there are two maps containing two categories: light and dark. We compute the proportion of cells that agree between the two maps and find that the agreement is 12/16 and the disagreement is 4/16. However, at a more sophisticated level the disagreement in terms of two components can be computed as (1) disagreement due to quantity, and (2) disagreement due to location. A disagreement due to quantity is defined as a disagreement between the maps in terms of the quantity of a category. For example, the proportion of cells in the dark category in the comparison map is 10/16 and in the reference map is 12/16; therefore, there is a disagreement due to quantity of 2/16. The example of disagreement of location could be when one replaces two cells from the comparison map with cells from the reference map to make them more consistent. Therefore, the difference in location is 2/16 [65].
In addition, fuzzy kappa was computed using the Map Comparison Kit (MCK) software [66]. The fuzzy kappa statistic is often used for the categorical maps with unequal legends (a different number of classes). The fuzzy kappa approach uses cells in the neighborhood to compute the equivalent of the traditional kappa statistic, but is instead based on fuzzy set theory [58,67,68]. One of three functions, which are exponential, linear, or constant, can be selected to determine the extent to which the neighboring cells influence the fuzzy representation [66]. Exponential decay (default setting) was chosen in the present research.
Finally, fraction correct and Cramer's V statistics were calculated. Fraction correct [68] is the number of equal cells divided by the total amount of cells. As an overall measure of similarity, the fraction correct method is considered a flawed estimate because it tends to consider maps with few or unevenly distributed categories as more similar than those with many, more equally distributed categories. Therefore, Cramer's V statistics were chosen to determine the strength of association between variables. Cramer's V statistics calculates the correlations in tables that have more than two rows and two columns. Values close to 0 denote a weak association between variables, while values close to 1 indicate a strong association.

LST-based typology consists of five main types of areas: (I) inland waters, (II) builtup areas, (III) agricultural/meadows and pastures, (IV) forested areas, and (V) mixed areas
(with a dominance of wetlands). The five types are divided into subtypes marked with small letters, for example, type IIIa ( Figure 9, Table 11), and described in terms of (1) land use/land cover (LULC), (2) lithology/moisture, (3) relief, and (4) normalized average land surface temperature (LST) (Tables 11 and 12).  * % of the area-descriptions of classes are generalized (percentage does not always add up to 100); only LULC, lithology/moisture, aspect/slope forms occupying more than 10% of a thermal subtype's area were included in the description. ** s.a.LST stands for standardized mean LST from all scenes/seasons used in the analysis.  Figure 4, the second column (topoclimatic type/subtype) refers to Figure 7: type I: inland waters; type II: built-up areas; type III: agricultural/meadows and pastures areas with subtypes a, b, and c; type IV: forested areas with subtypes a, b, c, and d; type V: mixed areas (with dominance of wetlands) with subtypes a and b.
Type I consists of inland waters such as watercourses, lakes, and reservoirs. Type II, characterized by high LST values (Table 12), occurs in urban and agricultural areas with a deep groundwater table (>1 m below ground surface). Sandy loam and sandy soil comprise the underlying geological materials. The largest percentage of this subtype underlies flat areas (~80%).
Type III consists of three subtypes. Subtype IIIa is entirely composed of agricultural flat areas, with underlying loam and sandy loam and a deep groundwater table (>1 m). The standardized average LST is moderate (Table 12). Subtype IIIb consists of mainly, but not only, agricultural flat areas. Underlying geological materials, however, are more heterogeneous compared to subtype IIIa, and composed of loam, sandy loam, and sandy soil. The groundwater table is deep and the standardized average LST is high. Subtype IIIc is composed of pastures and agricultural flat areas. The depth of the groundwater table varies depending on the underlying geological material, which is typically composed of peat/organic matter, sandy loam, and sandy soil. The standardized average LST is high.
Type IV is divided into four subtypes. Subtype IVa is made up of forested areas with a predominance of coniferous forests. Flat areas and slopes exposed to the north or west/east are dominant. Underlying soil consists of sand and sandy loam with a deep groundwater table (>1 m below the ground's surface). The standardized average LST is the lowest in forested areas of this subtype. Subtype IVb is composed of only forests, typically coniferous or mixed. Flat areas or slopes exposed to the west/east predominate. Soils consist of mainly of sandy loam with a deep groundwater table. The standardized average LST is low. Subtype IVc is composed of flat forested areas, where broadleaf forests prevail. The groundwater table varies with soil type, which is sandy loam and peat or organic matter. The standardized average LST is low but higher than in the IVa and IVb subtypes. Subtype IVd consists mainly of coniferous forests on predominately flat areas. Rare sloped areas exhibit a west or east aspect. The groundwater table is rather deep, with sandy loam in the subsurface. The standardized average LST is the highest in forested areas (Table 12).
Two subtypes constitute Type V. Subtype Va is covered by flat areas with coniferous forests, inland marshes, and pastures. Generally, the groundwater table is shallow, underlain by peat or organic matter and sandy soil. The standardized average LST temperature is moderate (Table 12). Subtype Vb consists mostly of pastures, complemented by coniferous forests and inland marshes. The area slopes less than in subtype Va, and exhibits a generally shallower groundwater table and higher humidity. The underlying soils have variable composition. The standardized, average LST temperature is moderate (Table 12). (Paszyński, 1980;Paszyński et al., 1999) The LST-based typology was compared with the Paszyński typology [28] using the cross-tabulation method (Tables 13 and 14). (1) From 41 to 59% of pixels in classes 11, 12, and 13 from Paszyński's map (urban and industrial areas) are contained in type II of the LST map (41%, 59%, and 53%, respectively), with the highest average standardized LST value of 1.50 (high agreement); (2) pixels in Paszyński's classes 5 and 6 (flat and concave areas) are distributed in different LST types/subtypes (low agreement); and (3) forested areas (classes 7, 8, 9, and 10 of Paszyński's map) have a relatively uniform, low, and balanced average standardized LST as compared to nonforested areas (Tables 11, 12, and 14). The KIA for the entire study area is 0.0015, which indicates a low overall agreement.

Comparison of the LST Based Typology with the Deductive Topoclimatic Classification
For the study area, Klocation and KlocationStrata (VALIDATE module) are both 0.0033, which indicates how well the grid cells are located on the landscape (Klocation) or within the strata/layer (KlocationStrata). The analysis of the agreement and disagreement for the study area shows that there is an 8% agreement due to chance, 4% agreement due to quantity, 1% agreement at the grid cell level, 43% disagreement at the grid cell level, and 49% disagreement due to quantity. For the kappa statistics, all values are >0%, where 0% indicates that the level of agreement is equal to the agreement due to chance and 100% indicates perfect agreement [65]. The overall agreement between the two maps is weak, but by analyzing the smaller areas and particular classes, one can see that its distribution includes areas with a very high agreement between the two maps.
In addition, two fuzzy kappa maps were constructed. The first ( Figure 10) is based on the 0/1 contingency table (Table 13), where 1 is assigned to the consistent classes. The fuzzy kappa statistic is 0.011 and the fraction correct is 0.619. The second map (Figure 11) was created from the proportional cross tabulation table (Table 14). Here, the fuzzy kappa statistic is 0.063, whereas the fraction correct is 0.214. Overall, both statistics present low values, which means that the agreement between the two classifications is weak. However, the fuzzy kappa maps provide important information on the spatial distribution of kappa statistics for Wolin Island.   Analyzing the MCK fuzzy kappa maps ( Figures 10 and 11), the areas with high agreement are localized. For example, water bodies and urban areas (distinguished by green/yellow color in Figure 10) have the highest fuzzy kappa values (0.9-1). In contrast, lowlands and wetlands have the lowest kappa values (0-0.2). Terminal moraines covered by forest areas are characterized by moderate fuzzy kappa values (0.3-0.4). These findings are consistent with the cross-tabulation results (Table 14). For instance, water bodies and urban areas yielded by cross-tabulation analysis have the highest fuzzy kappa values on both maps. They are clearly distinguished by a green/yellow color in Figure 10, and shown even more clearly in Figure 11. Cramer's V value, which determines the strength of the relationship between two maps, equals 0.38. This indicates that there is weak correlation between the distributions of classes on both maps. However, specified pairs of classes show similarities (e.g., urban areas, water), as indicated by the MCK maps (Figures 10 and 11).
Arable lands (211, 222, 231, 242, and 243 in Figure 2) show a value of 1 on the 0/1 cross tabulation map (Table 13, Figure 10), which means that the classes of the LST map and Paszyński's map overlap. However, one category in Paszyński's classification (Figure 8class 5: flat forms) encloses seven thermal classes on the LST map ( Figure 9).
Since the LSTs differ significantly in LST classes, even in flat areas, it can be hypothesized that Paszyński's classification should be more detailed for flat areas. In fact, Ryszkowki and Kędziora [4] proved using in situ measurements that the values of particular components of heat balance in flat areas vary significantly.
The LST-based topoclimate typology was compared with Paszyński's classification. It was found that (1) water bodies are characterized by the lowest temperatures, (2) unforested areas with specific exposures belong to the same classes on the LST map, (3) forested areas are characterized by more balanced, relatively low LST values on LST maps, (4) pixels of flat and concave areas (from Paszyński) are distributed in several thermal classes without one dominant class, and (5) urban and industrial classes from Paszynski's classification have the highest LST.
Traditional methods (Paszyński et al., in this work) show topoclimatic differentiation in a topographically diverse area, as they are based mainly on relief characteristics. We focused on flat areas, as such classifications are often demanded by local authorities and farmer associations to optimize agricultural production. The advantage of the technique presented is that it shows good sensitivity to changes in topoclimatic conditions (demonstrated by changes in LST distribution) even in flat agricultural areas with little topographic variability. However, it also gives reliable results on moraine hills and cliffs occurring on Wolin Island (Areas 1 and 3, Figure 1).
The constraint of the presented method is the fact that research based on satellite data are always limited in terms of temporal or spatial resolution. The method presented can be further developed and modified by the use of the other data, e.g., atmospheric properties from numerical weather prediction models (European Centre for Medium-Range Weather Forecasts, or ECMWF, for instance).