Next Article in Journal
Mapping Threats of Spring Frost Damage to Tea Plants Using Satellite-Based Minimum Temperature Estimation in China
Previous Article in Journal
Gradient Boosting Machine and Object-Based CNN for Land Cover Classification
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Research Centre in Wrocław, Institute of Geological Sciences, Polish Academy of Sciences, Podwale 75, 50-449 Wrocław, Poland
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(14), 2712; https://doi.org/10.3390/rs13142712
Submission received: 25 May 2021 / Revised: 6 July 2021 / Accepted: 7 July 2021 / Published: 9 July 2021

Abstract

:
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 topography. The technique will hopefully prove to be a convenient and relatively fast tool that can improve the topoclimatic classification of other regions. It could be applied by local authorities and farmer associations for optimizing agricultural production.

1. Introduction

Local climates are controlled by global, regional, and local factors. For example, planetary waves operate in large to macro scales (>103 km) and high-pressure systems or trajectories of cyclones contribute to mesoscale processes (101–103 km). The most common variation in climatic patterns, however, is at topoclimatic scales between 101 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).
Topoclimate studies are mostly based on in situ measurements [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]. In this work, a new, alternative method for topoclimate determination over a large area (the extent >10 km) based on land surface temperature (LST) distribution from thermal remote sensing (TRS) data is presented. At the local scale, LST is one of the factors influencing and describing topoclimatic diversity. LST has a direct impact on air temperature and is one of the main physical characteristics reflecting the processes on the earth’s surface [19].
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.
The results of the LST-based topoclimate mapping of Wolin Island supplement the existing body of knowledge on climate and thermal remote sensing applications, including (1) vegetation studies [29]; (2) soil moisture assessments [30]; (3) topoclimatic research [5,31]; (4) thermal anomalies [32]; (5) permafrost detection [33]; and (6) the detection of urban heat islands [21,34,35,36]. These works all point to the importance of remotely sensed LSTs in environmental studies. The main advantage of satellite images is that they provide simultaneous data for relatively large areas.

2. Study Area

Wolin, the largest island (265 km2) 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. Non-irrigated 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).

3. 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).

3.1. 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).
Table 1. Circulation types during the Landsat ETM+ images acquisition times. Circulation types are according to [41].
Table 1. Circulation types during the Landsat ETM+ images acquisition times. Circulation types are according to [41].
DateAcquisition Time (GMT)Circulation TypeGroupDuration (Days)
11.07.1999 *09:55Fennoscandian high, cyclonicmeridional3
13.09.199909:55Central European highmixed6
08.04.200009:54British Islands highmeridional2
14.08.200009:53Central European ridgemixed9
13.05.200109:52Norwegian Sea–Fennoscandian high, anticyclonicmeridional6
20.08.200209:50Fennoscandian high, anticyclonicmeridional6
16.03.200309:51British Islands highmeridional8
* Landsat images IDs: 11.07.1999—LE71930221999192EDC00, 13.09.1999—LE71930221999256EDC00, 08.04.2000—LE71930222000099EDC00, 14.08.2000—LE71930222000227EDC00, 13.05.2001—LE71930222001133EDC00, 20.08.2002—LE71930222002232EDC00, 16.03.2003—LE71930222003075EDC00.
Figure 5. Seven LST (land surface temperature) images of Wolin Island calculated using the LST program [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.
Figure 5. Seven LST (land surface temperature) images of Wolin Island calculated using the LST program [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.
Remotesensing 13 02712 g005

3.2. 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:
T s = a 6 ( 1 C 6 D 6 ) + [ b 6 ( 1 C 6 D 6 ) + C 6 + D 6 ] T s e n s o r D 6 T a C 6 ,
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:
C6 = ετ6,
D6 = (1 − τ6)[1 + (1 − ε)τ6],
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).
Ta is calculated from one of linear Equations (4)–(7) referring to the four standard atmospheres:
Ta = 1.9769 + 0.91715T0 (for USA 1976),
Ta = 17.9769 + 0.91715T0 (for tropical),
Ta = 16.0110 + 0.92621T0 (for mid-latitude summer),
Ta = 19.2704 + 0.91118T0 (for mid-latitude winter),
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:
w = w z R W z ,
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
g z = 219.5 × e T
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]:
f = e e m a x 100 % .
Therefore,
e = f 100 % e m a x ,
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:
NDVI = (band4 − band3)/(band4 + band3).
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 (Figure 1 and Figure 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.

3.3. 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.

3.4. Factors Influencing the LST Spatial Variability

Thermal classes were described in terms of (1) LULC, (2) lithology/moisture, (3) relief, and (4) average standardized LST.
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.
Relief characteristics were calculated using Digital Terrain Elevation Data Level 2 (DTED2) [60]. Three classes of terrain slopes were defined in accordance with [61]: (1) flat areas (<2°), (2) low slopes (2–5°), and (3) high slopes (>5°). Aspects were categorized into four classes: (1) south (135–225°), (2) east and west (45–135° and 225–315°), (3) north (0–45° and 315–360°), and (4) flat areas without aspect. The slope and aspect classes were then merged and recoded into seven combined classes (Table 9 and Figure 3).

3.5. 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,
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 LST-based 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).

3.6. 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 grid-cell 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.

4. Results

LST-based typology consists of five main types of areas: (I) inland waters, (II) built-up 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) (Table 11 and Table 12).
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).

5. Discussion

Comparison of the LST Based Typology with the Deductive Topoclimatic Classification (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 (Table 13 and Table 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 (Table 11, Table 12, and Table 14). The KIA for the entire study area is 0.0015, which indicates a low overall agreement.
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 (Figure 10 and Figure 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 (Figure 10 and Figure 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 8—class 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.

6. Conclusions

The proposed LST-based typology includes five basic types of topoclimate: (I) inland waters, (II) built-up areas, (III) agricultural/meadows and pastures, (IV) forested areas, and (V) mixed areas (with a dominance of wetlands). The five types are subdivided according to (1) the land use/land cover (LULC), (2) lithology/moisture, (3) relief, and (4) normalized average land surface temperature (LST).
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).

Author Contributions

Conceptualization, M.C.; methodology, M.C., J.C.; software, M.C.; validation, M.C.; formal analysis, M.C.; investigation, M.C.; resources, M.C.; data curation, M.C.; writing—original draft preparation, M.C.; writing—review and editing, J.C.; visualization, M.C.; project administration, M.C.; funding acquisition, M.C., J.C. Both authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Ministry of Science, Innovation and Higher Education, grant number N N306 140538. Jakub Ciążela is funded from the National Science Centre Poland, program OPUS 19, grant no. 2020/37/B/ST10/01420.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

This research was supported by a grant from Ministry of Science, Innovation and Higher Education entitled: “Effect of thermal properties of Earth surface on water circulation and biogeochemical processes in the last glacial period” (no. N N306 140538). Jakub Ciążela is funded from the National Science Centre Poland, program OPUS 19, grant no. 2020/37/B/ST10/01420. M.Ciazela thanks Alfred Stach for his help in designing this study.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Slavich, E.; Warton, D.I.; Ashcroft, M.B.; Gollan, J.R.; Ramp, D. Topoclimate versus macroclimate: How does climate mapping methodology affect species distribution models and climate change projections? Divers. Distrib. 2014, 20, 952–963. [Google Scholar] [CrossRef]
  2. Lembrechts, J.J.; Lenoir, J.; Roth, N.; Hattab, T.; Milbau, A.; Haider, S.; Pellissier, L.; Pauchard, A.; Backes, A.R.; Dimarco, R.D.; et al. Comparing temperature data sources for use in species distribution models: From in-situ logging to remote sensing. Glob. Ecol. Biogeogr. 2019, 28, 1578–1596. [Google Scholar] [CrossRef]
  3. Kluge, M. Methods of Topoclimatic Studies. In Dokumentacja Geograficzna (Geographic Documentation); Wydawnictwo Polskiej Akademii Nauk: Warszawa, Poland, 1980. [Google Scholar]
  4. Ryszkowski, L.; Kędziora, A. Impact of agricultural landscape structure on energy flow and water cycling. Landsc. Ecol. 1987, 1, 85–94. [Google Scholar] [CrossRef] [Green Version]
  5. Wojkowski, B.; Skowera, J. The structure of the radiation balance dependent on synoptic situations. Pamiętnik Puławski 2007, 144, 169–178. [Google Scholar]
  6. Cyrulewski, K.; Suwała, K. Diversity of air temperature and humidity in the chosen part of Różany Potok river catchment. Badania Fizjograficzne 2008, 59, 43–52. [Google Scholar] [CrossRef]
  7. Dobrowski, S.Z.; Abatzoglou, J.T.; Greenberg, J.A.; Schladow, S.G. How much influence does landscape-scale physiography have on air temperature in a mountain environment? Agric. For. Meteorol. 2009, 149, 1751–1758. [Google Scholar] [CrossRef]
  8. Holden, Z.A.; Crimmins, M.A.; Cushman, S.A.; Littell, J.S. Empirical modeling of spatial and temporal variation in warm season nocturnal air temperatures in two North Idaho mountain ranges, USA. Agric. For. Meteorol. 2011, 151, 261–269. [Google Scholar] [CrossRef]
  9. Puşcoi, G.; Murărescu, B.; Pehoiu, O. Aspects of Topoclimate Regionalization in Targoviste City (Romania). In Water Resources and Wetlands; Valahia University of Targoviste: Targoviste, Romania, 2012. [Google Scholar]
  10. Convey, P.; Coulson, S.J.; Worland, M.R.; Sjöblom, A. The importance of understanding annual and shorter-term temperature patterns and variation in the surface levels of polar soils for terrestrial biota. Polar Biol. 2018, 41, 1587–1605. [Google Scholar] [CrossRef] [Green Version]
  11. Grzybowski, J. Attempt to Distinguish Topoclimates Types in Poland. In Problemy Współczesnej Topoklimatologii (Problems of Contemporary Topoclimatology); IGiPZ PAN: Warszawa, Poland, 1990; pp. 34–40. [Google Scholar]
  12. Tamulewicz, J. Radiation Balance of the Agricultural Landscape on the Example of the Vicinity of Turwia. In Problemy Współczesnej Topoklimatologii (Problems of Contemporary Topoclimatology); IGiPZ PAN: Warszawa, Poland, 1990; pp. 69–78. [Google Scholar]
  13. Błażejczyk, K.; Krawczyk, B.; Skoczek, J. Badania topoklimatyczne i mikroklimatyczne w różnych strefach klimatycznych = Topoclimatic and microclimatic investigations in different climatic zones. In Zeszyty Instytutu Geografii i Przestrzennego Zagospodarowania PAN; IGiPZ PAN: Warszawa, Poland, 1992; Volume 5. [Google Scholar]
  14. Ryszkowski, A.; Kędziora, L. Modification of the effects of global climate change by plant cover structure in an agricultural landscape. Geogr. Pol. 1995, 65, 5–35. [Google Scholar]
  15. Richards, K.; Baumgarten, M. Towards Topoclimate Maps of Frost and Frost Risk for Southland, New Zealand. In Proceedings of the 15th Annual Colloquium of the Spatial Information Research Centre (SIRC 2003: Land, Place and Space), Dunedin, New Zealand, 1–2 December 2003. [Google Scholar]
  16. Eulenstein, F.; Urbaniak, M.; Chojnicki, B.H.; Olejnik, J. Influence of plant cover on the share of the soil heat flux in the heat balance of the active surface. Int. Agrophys. 2005, 19, 31–36. Available online: http://www.international-agrophysics.org/Influence-of-plant-cover-on-the-share-of-the-soil-heat-flux-in-the-heat-balance-of,106662,0,2.html (accessed on 26 June 2021).
  17. Vysoudil, M.; Navrátil, L. Topoclimatological research in Údolí Bystřice River Nature Park (Czech Republic): Functional meteorological network. Geographica 2006, 39, 111–139. [Google Scholar]
  18. Kossowski, J.; Łykowski, B. Daily sums of solar radiation during summer period in Felin near Lublin and their relationship with sunshine duration and cloudiness. Sci. Rev. Eng. Environ. Sci. 2007, 2, 74–84. [Google Scholar]
  19. Zhang, J.; Wang, Y.; Li, Y. A C++ program for retrieving land surface temperature from the data of Landsat TM/ETM+ band6. Comput. Geosci. 2006, 32, 1796–1805. [Google Scholar] [CrossRef]
  20. Balew, A.; Korme, T. Monitoring land surface temperature in Bahir Dar city and its surrounding using Landsat images. Egypt J. Remote Sens. Space Sci. 2020, 23, 371–386. [Google Scholar] [CrossRef]
  21. Cai, G.; Du, M.; Xue, Y. Monitoring of urban heat island effect in Beijing combining ASTER and TM data. Int. J. Remote Sens. 2011, 32, 1213–1232. [Google Scholar] [CrossRef]
  22. Clinton, N.; Gong, P. MODIS detected surface urban heat islands and sinks: Global locations and controls. Remote Sens. Environ. 2013, 134, 294–304. [Google Scholar] [CrossRef]
  23. Chen, X.; Zhang, Y. Impacts of urban surface characteristics on spatiotemporal pattern of land surface temperature in Kunming of China. Sustain. Cities Soc. 2017, 32, 87–99. [Google Scholar] [CrossRef] [Green Version]
  24. Sobrino, J.A.; Jiménez-Muñoz, J.C.; Paolini, L. Land surface temperature retrieval from LANDSAT TM 5. Remote Sens. Environ. 2004, 90, 434–440. [Google Scholar] [CrossRef]
  25. Sobrino, J.A.; Oltra-Carrió, R.; Sòria, G.; Jiménez-Muñoz, J.C.; Franch, B.; Hidalgo, V.; Mattar, C.; Julien, Y.; Cuenca, J.; Romaguera, M.; et al. Evaluation of the surface urban heat island effect in the city of Madrid by thermal remote sensing. Int. J. Remote Sens. 2013, 34, 3177–3192. [Google Scholar] [CrossRef]
  26. Stathopoulou, M.; Cartalis, C. Daytime urban heat islands from Landsat ETM+ and Corine land cover data: An application to major cities in Greece. Sol. Energy 2007, 81, 358–368. [Google Scholar] [CrossRef]
  27. Yuan, F.; Bauer, M.E. Comparison of impervious surface area and normalized difference vegetation index as indicators of surface urban heat island effects in Landsat imagery. Remote Sens. Environ. 2007, 106, 375–386. [Google Scholar] [CrossRef]
  28. Paszyński, J.; Miara, K.; Skoczek, J. Topoclimatological mapping as based on the energy exchange at the interface Earth-atmosphere. In Dokumentacja Geograficzna (Geographical Documentation); IGiPZ PAN: Warszawa, Poland, 1999; Volume 14. [Google Scholar]
  29. Berni, J.A.J.; Zarco-Tejada, P.J.; Suárez, L.; Fereres, E. Thermal and narrowband multispectral remote sensing for vegetation monitoring from an unmanned aerial vehicle. IEEE Trans. Geosci. Remote Sens. 2009, 47, 722–738. [Google Scholar] [CrossRef] [Green Version]
  30. Hermanowska, B. Numeryczne Modelowanie Inercji Termalnej Gruntu dla Teledetekcyjnego Określania Jego Wilgotności; Akademia Górniczo-Hutnicza im. Stanisława Staszica: Kraków, Poland, 1997. [Google Scholar]
  31. Kubiak, M.; Dzieszko, P. Using Thermal Remote Sensing in Environmental Studies. Trans. GIS 2012, 16, 715–732. [Google Scholar] [CrossRef]
  32. Kuenzer, C.; Hecker, C.; Zhang, J.; Wessling, S.; Wagner, W. The potential of multidiurnal MODIS thermal band data for coal fire detection. Int. J. Remote Sens. 2008, 29, 923–944. [Google Scholar] [CrossRef]
  33. Kędzia, S.; Mościcki, J. Geomorphology of the Carpatho-Balcan Region. In The multiannual permafrost occurence in the High Tatra Mountains Kozia Valley, Proceedings of the Carpatho-Balcan Conference, Baile Herculane, Orsova, Drobeta Turnu Severin, Romania, 11–17 October 1998; Case study; IGiPZ PAN: Warszawa, Poland, 2000; pp. 100–106. [Google Scholar]
  34. Osińska-Skotak, K.; Madany, A. Use of LANDSAT TM satellite data to determine the characteristics of the Warsaw urban heat island. Prace Naukowe Politechniki Warszawskiej. Inżynieria Środowiska (Scientific Papers of the Warsaw University of Technology-Environmental Engineering) 1998, 26, 6–33. (In Polish). Available online: https://www.academia.edu/21232622/Wykorzystanie_danych_satelitarnych_LANDSAT_TM_do_okre%C5%9Blenia_warszawskiej_wyspy_ciep%C5%82a (accessed on 9 July 2021).
  35. Voogt, J.A.; Oke, T.R. Thermal remote sensing of urban climates. Remote Sens. Environ. 2003, 86, 370–384. [Google Scholar] [CrossRef]
  36. Amarsaikhan, D.; Ganzori, M.; Tae-Heon, M. Investigation of Urban Tem-Perature Changes Using Multitemporal Thermal Infrared Images. In Proceedings of the 26th Asian Conference on Remote Sensing and 2nd Asian Space Conference (ACRS2005), Hanoi, Vietnam, 7–11 November 2005. [Google Scholar]
  37. Kostrzewski, A. Studia z Geografii Fizycznej i Ekonomicznej Wyspy Wolin; SKNG UAM: Poznan, Poland, 1978. [Google Scholar]
  38. Kottek, M.; Grieser, J.; Beck, C.; Rudolf, B.; Rubel, F. World map of the Köppen-Geiger climate classification updated. Meteorologische Zeitschrif 2006, 15, 259–263. [Google Scholar] [CrossRef]
  39. Rubel, F.; Kottek, M. Observed and projected climate shifts 1901-2100 depicted by world maps of the Köppen-Geiger climate classification. Meteorologische Zeitschrift 2010, 19, 135–141. [Google Scholar] [CrossRef] [Green Version]
  40. Rubel, F.; Brugger, K.; Haslinger, K.; Auer, I. The climate of the European Alps: Shift of very high resolution Köppen-Geiger climate zones 1800–2100. Meteorologische Zeitschrift 2017, 26, 115–125. [Google Scholar] [CrossRef]
  41. Hess, P.; Brezowsky, H. Katalog der Grosswetterlagen Europas. Neubearbaitete und ergantze Auflage. Berichte des Deutschen Wetterdienstes; Potsdam Institute for Climate Impact Research: Potsdam, Germany, 1969; Volume 113, p. 68. [Google Scholar]
  42. Qin, Z.; Karnieli, A.; Berliner, P. A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region. Int. J. Remote Sens. 2001, 22, 3719–3746. [Google Scholar] [CrossRef]
  43. Snyder, W.C.; Wan, Z.; Zhang, Y.; Feng, Y.-Z. Classification-based emissivity for land surface temperature measurement from space. Int. J. Remote Sens. 1998, 19, 2753–2774. [Google Scholar] [CrossRef]
  44. Paszyński, J. Methods of topoclimatic mapping. Dok. Geogr. Geogr. Doc. 1980, 3, 13–28. [Google Scholar]
  45. Richards, J.A. Remote Sensing Digital Image Analysis: An Introduction; Springer: Berlin/Heidelberg, Germany, 2013. [Google Scholar]
  46. ScienceOpen Inc. STATISTICA (data analysis software system), version 8.0; ScienceOpen Inc.: Burlington, MA, USA, 2008. [Google Scholar]
  47. Ball, G.H.; Hall, D.J. Isodata, A Novel Method of Data Analysis and Pattern Classification. 1965. Available online: http://www.dtic.mil/cgi-bin/GetTRDoc?Location=U2&doc=GetTRDoc.pdf&AD=AD0699616 (accessed on 7 April 2021).
  48. Hartigan, J.A.; Wong, M.A. Algorithm AS 136: A K-Means Clustering Algorithm. Appl. Stat. 1979, 28, 100. [Google Scholar] [CrossRef]
  49. Kubiak, M.; Stach, A. The use of motor-glider in topoclimatic studies. Int. J. Appl. Earth Obs. Geoinf. 2014, 32, 86–198. [Google Scholar] [CrossRef]
  50. Matkowska, Z.; Ruszała, M.; Wdowiak, M. Explanatory Booklet to the Detailed Geological Map of Poland 1:50,000 Sheet Międzyzdroje; PGI: Warsaw, Poland, 1977. Available online: http://bazadata.pgi.gov.pl/data/smgp/arkusze_txt/smgp0113.pdf (accessed on 9 July 2021).
  51. Ruszała, M.; Wdowiak, M. Detailed Geological Map of Poland 1:50,000 sheet Międzyzdroje no. 113/1977. In Wydaw. Geol. Warszawa; 1977. Available online: http://bazadata.pgi.gov.pl/data/smgp/arkusze_skany/smgp0113.jpg (accessed on 9 July 2021).
  52. Ruszała, M. Detailed geological map of Poland 1:50,000, sheet Wolin no. 114/1978. In Wydaw. Geol. Warszawa; 1978. Available online: http://bazadata.pgi.gov.pl/data/smgp/arkusze_skany/smgp0114.jpg (accessed on 9 July 2021).
  53. Ruszała, M. Detailed Geological Map of Poland 1:50,000 sheet Racimierz no. 152/198. In Wydaw. Geol. Warszawa; 1981. Available online: http://bazadata.pgi.gov.pl/data/smgp/arkusze_skany/smgp0152.jpg (accessed on 9 July 2021).
  54. Ruszała, M. Explanatory booklet to the Detailed Geological Map of Poland 1:50,000 sheet Racimierz. In Wydaw. Geol. Warszawa; 1981. Available online: https://bazadata.pgi.gov.pl/data/smgp/arkusze_txt/smgp0152.pdf (accessed on 9 July 2021).
  55. Ruszała, M.; Dobracka, E.; Piotrowski, A. Explanatory booklet to the Detailed Geological Map of Poland 1:50,000 sheet Międzywodzie. In Wydaw. Geol. Warszawa; 1979. Available online: https://bazadata.pgi.gov.pl/data/smgp/arkusze_txt/smgp0114.pdf (accessed on 9 July 2021).
  56. Chesworth, W. Encyclopedia of Soil Science; Springer: New York, NY, USA, 2008. [Google Scholar]
  57. Arya, S.P. Introduction to Micrometeorology. In International Geophysics Series; Academic Press: Cambridge, MA, USA, 2001. [Google Scholar]
  58. Rose, D.A.; Lal, R.; Shukla, K.M. Principles of Soil Physics; CRC Press: Boca Raton, FL, USA, 2004; ISBN 0-8247-5324-0. [Google Scholar]
  59. GUGiK (Główny Urząd Geodezji i Kartografii). Technical Guidelines GIS—Hydrographic Map of Poland in the Scale 1:50 000; GUGiK (Główny Urząd Geodezji i Kartografii): Warsaw, Poland, 2005. [Google Scholar]
  60. National Imagery and Mapping Agency. Performance Specification Digital Terrain Elevation Data (DTED); National Imagery and Mapping Agency: Springfield, VA, USA, 2000.
  61. Klimaszewski, M. Geomorphology; Polish Scientific Publishers: Warsaw, Poland, 1978. [Google Scholar]
  62. Jasiewicz, J.; Stepinski, T.F. Geomorphons-a pattern recognition approach to classification and mapping of landforms. Geomorphology 2013, 182, 147–156. [Google Scholar] [CrossRef]
  63. Neteler, M.; Mitasova, H. Open Source GIS: A GRASS GIS Approach; Springer US: New York, NY, USA, 2008. [Google Scholar]
  64. Cohen, J. Weighted kappa: Nominal scale agreement with provision for scaled disagreement or partial credit. Psychol. Bull. 1968, 70, 213–220. [Google Scholar] [CrossRef]
  65. Pontius, R.G. Statistical methods to partition effects of quantity and location during comparison of categorical maps at multiple resolutions. Photogramm. Eng. Remote Sens. 2002, 68, 1041–1050. [Google Scholar]
  66. Visser, H.; de Nijs, T. The map comparison kit. Environ. Model. Softw. 2006, 21, 346–358. [Google Scholar] [CrossRef]
  67. Hagen, A. Fuzzy set approach to assessing similarity of categorical maps. Int. J. Geogr. Inf. Sci. 2010, 17, 235–249. [Google Scholar] [CrossRef] [Green Version]
  68. Vliet, J. MCK User Manual; Research Institute for Knowledge Systems (RIKS BV): Maastricht, The Netherlands, 2010. [Google Scholar]
Figure 1. Digital terrain elevation data level 2 of Wolin Island with the most characteristic landscape features (marked with numbers): (1) moraine hills covered with forest, mainly beech (Pomeranian beech); (2) swampy and lowland areas, such as a trench peninsula in the south of the island; (3) the highest coastal cliffs in Poland with a height of about 80 m extending over 15 km; (4) beaches; (5) the reverse delta of Świna River–Świna Gate; (6) and a lake district consisting of nine lakes and Lewińska stream, which drains water from the lake district to the Dźiwna River.
Figure 1. Digital terrain elevation data level 2 of Wolin Island with the most characteristic landscape features (marked with numbers): (1) moraine hills covered with forest, mainly beech (Pomeranian beech); (2) swampy and lowland areas, such as a trench peninsula in the south of the island; (3) the highest coastal cliffs in Poland with a height of about 80 m extending over 15 km; (4) beaches; (5) the reverse delta of Świna River–Świna Gate; (6) and a lake district consisting of nine lakes and Lewińska stream, which drains water from the lake district to the Dźiwna River.
Remotesensing 13 02712 g001
Figure 2. The CORINE Land Cover 2006 (CLC2006) map (http://www.eea.europa.eu/publications/COR0-part1, accessed on 26 June 2021) of land use/land cover forms.
Figure 2. The CORINE Land Cover 2006 (CLC2006) map (http://www.eea.europa.eu/publications/COR0-part1, accessed on 26 June 2021) of land use/land cover forms.
Remotesensing 13 02712 g002
Figure 3. Aspect/slope classes of Wolin Island.
Figure 3. Aspect/slope classes of Wolin Island.
Remotesensing 13 02712 g003
Figure 4. Research scheme.
Figure 4. Research scheme.
Remotesensing 13 02712 g004
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,48].
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,48].
Remotesensing 13 02712 g006
Figure 7. Lithology/moisture classes of Wolin Island.
Figure 7. Lithology/moisture classes of Wolin Island.
Remotesensing 13 02712 g007
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.
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.
Remotesensing 13 02712 g008
Figure 9. Proposed thermal typology of Wolin Island based on the spatial distribution of LST (derived from Landsat 7 thermal images). 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.
Figure 9. Proposed thermal typology of Wolin Island based on the spatial distribution of LST (derived from Landsat 7 thermal images). 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.
Remotesensing 13 02712 g009
Figure 10. Fuzzy kappa map based on the 0/1 contingency table (Table 13). A value of 1 (green) means total similarity and a value of 0 (red) means total dissimilarity. Values in yellow (around 0.5) indicate some similarity between the Paszyński (Figure 8) classification and LST (Figure 9) maps.
Figure 10. Fuzzy kappa map based on the 0/1 contingency table (Table 13). A value of 1 (green) means total similarity and a value of 0 (red) means total dissimilarity. Values in yellow (around 0.5) indicate some similarity between the Paszyński (Figure 8) classification and LST (Figure 9) maps.
Remotesensing 13 02712 g010
Figure 11. Fuzzy kappa map based on a proportional cross-tabulation table (Table 14). Values range from 0 to 0.79, where 1 (green) means total similarity and 0 (red) means total dissimilarity. Values in orange (around 0.5) indicate some similarity between the Paszyński (Figure 8) classification and LST (Figure 9) maps.
Figure 11. Fuzzy kappa map based on a proportional cross-tabulation table (Table 14). Values range from 0 to 0.79, where 1 (green) means total similarity and 0 (red) means total dissimilarity. Values in orange (around 0.5) indicate some similarity between the Paszyński (Figure 8) classification and LST (Figure 9) maps.
Remotesensing 13 02712 g011
Table 2. Estimation of atmospheric transmittance for Landsat ETM+ channel 6 [42]. Band 6L stands for low gain and 6H stands for high gain.
Table 2. Estimation of atmospheric transmittance for Landsat ETM+ channel 6 [42]. Band 6L stands for low gain and 6H stands for high gain.
Band 6 ETM+Water Vapor (w) (g cm−2)Transmittance Estimation Equation
Band 6L0.4–1.6τ6 = 0.974–0.08007w
1.6–3.0τ6 = 1.031–0.1153w
Band 6H0.4–1.6τ6 = 0.982–0.09611w
1.6–3.0τ6 = 1.054–0.141w
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].
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].
Altitude (km)USA 1976TropicalMid-Latitude SummerMid-Latitude WinterAverage Rw(z)
00.4020580.4250430.4384460.4001240.416418
10.2562340.2610320.2621000.2542100.258394
20.1583230.1684000.1489430.1618730.159385
30.0874950.0759990.0744710.0955280.083373
40.0474970.0318780.0383640.0465100.041062
50.0245120.0193810.0179250.0237110.021382
60.0128460.0097710.0097360.0115140.010967
70.0062500.0047820.0052230.0040920.005087
80.0031320.0022570.0026110.0014710.002368
90.0010490.0009540.0013150.0005870.000976
100.0003580.0003490.0006160.0002380.000390
110.0001420.0001040.0001850.0000600.000123
120.0000550.0000320.0000440.0000260.000039
130.0000230.0000080.0000090.0000160.000014
140.0000090.0000040.0000040.0000110.000007
150.0000060.0000020.0000020.0000080.000004
Table 4. Correlation coefficients of seven LST images calculated with the Qin et al. algorithm [42].
Table 4. Correlation coefficients of seven LST images calculated with the Qin et al. algorithm [42].
Date11.07.199913.09.199908.04.200014.08.200013.05.200120.08.200216.03.2003
11.07.19991.000.680.690.790.730.690.65
13.09.1999 1.000.650.830.560.880.36
08.04.2000 1.000.670.780.660.65
14.08.2000 1.000.670.820.55
13.05.2001 1.000.550.71
20.08.2002 1.000.37
16.03.2003 1.00
Table 5. Single-factor analysis of the variance in the classes of the LST classification map based on the isodata algorithm (ISODATA11).
Table 5. Single-factor analysis of the variance in the classes of the LST classification map based on the isodata algorithm (ISODATA11).
EffectTest *ValueFEffect dfError dfp
InterceptWilks0.9034716.673089210.000001
ISODATA11Wilks0.00537267.37018013100.000001
* Wilks—Wilks’s lambda test; F—F-statistic; df—degrees of freedom; p—probability.
Table 6. Post hoc tests results. The number of combinations (from 7 selected Landsat images and comparing each image with every other one) was insignificantly different (similar) from the total number of combinations. The numbers in the map names stand for the number of classes. The calculations were carried out using the IDRISI software, except for the ISODATA8_ENVI map.
Table 6. Post hoc tests results. The number of combinations (from 7 selected Landsat images and comparing each image with every other one) was insignificantly different (similar) from the total number of combinations. The numbers in the map names stand for the number of classes. The calculations were carried out using the IDRISI software, except for the ISODATA8_ENVI map.
MapTotal Number of CombinationsCombinations Not Significantly Different%
ISOCLUST734310.3
ISOCLUST15157590.6
ISODATA1184740.5
ISODATA8_ENVI44800
CLUSTER1008111.1
Table 7. Discriminant analysis summary table based on the LST classifications maps calculated with the Qin et al. algorithm. The numbers in the map names stand for the number of classes, ENVI refers to ENVI software, and no ENVI indicates IDRISI software.
Table 7. Discriminant analysis summary table based on the LST classifications maps calculated with the Qin et al. algorithm. The numbers in the map names stand for the number of classes, ENVI refers to ENVI software, and no ENVI indicates IDRISI software.
Classification MethodMapNo of ClassesWilks’ LambdaPercent Correct
isoclustISOCLUST87 *0.048183%
ISOCLUST1615 *0.010379%
isodataISODATA11110.005494%
ISODATA8_ENVI80.011293%
clusterCLUSTER12120.019582%
* One of the classes is background.
Table 8. Lithology/moisture classes for Wolin Island.
Table 8. Lithology/moisture classes for Wolin Island.
ClassDescriptionArea (km2)% of Total Area
1Small water bodies3.31.1
2Peat/organic matter—wet61.721.6
3Sandy soil—wet24.88.7
4Sandy loam—wet15.85.5
5Loam/clay loam/silt loam—wet5.31.9
6Clay—wet0.30.1
7Peat/organic matter—dry4.61.6
8Sandy soil—dry42.915.0
9Sandy loam—dry99.634.8
10Loam/clay loam/silt loam—dry20.87.3
11Lakes6.82.4
Total:285.9100
Table 9. Aspect/slope classes.
Table 9. Aspect/slope classes.
ClassDescriptionArea (km2)% of Wolin Island’s Area
1No aspect areas/flat areas206.274.0
2South/low slope areas (2–5°)12.54.5
3East–west/low slope areas (2–5°)21.37.6
4North/low slope areas (2–5°)11.74.2
5South/high slope areas (>5°)7.92.8
6East–west/high slopes areas (>5°)11.94.3
7North/high slope areas (>5°)7.22.6
Table 10. Topoclimatic classification based on modified Paszyński’s theoretical assumptions [28,44].
Table 10. Topoclimatic classification based on modified Paszyński’s theoretical assumptions [28,44].
IDClassDescriptionNew Code
1.Convex forms
11.1.South aspect and slope > 5°2
21.2.All areas exclusive north and south aspects and slopes > 5°3
Areas with north and south aspects and slopes < 5°
31.3.North aspect and slope > 5°4
2.Flat forms, slopes < 5°5
3.Concave forms6
4.Forested areas
104.1.South exposition, slope > 5°7
114.2.Forested flat forms8
All areas with exclusive north or south aspects, slope > 5°9
124.3.North aspect, slope > 5°10
5.Urban and industrial areas
135.1.With convex forms11
145.2.With plains12
155.3.With valleys13
6.Water1
Table 11. Description of the proposed thermal types and subtypes on Wolin Island based on the spatial distribution of LST derived from Landsat 7 images.
Table 11. Description of the proposed thermal types and subtypes on Wolin Island based on the spatial distribution of LST derived from Landsat 7 images.
Topoclimate TypeLULCLithology/MoistureReliefs.a.LST (°C) **
NameCLC Code% *LithologyGroundwater%DescriptionInclination%
Type I—Inland waterswater courses511 −1.91
water bodies512
Type II -built-up areasartificial surfaces112, 12330sandy loam>1 m50flat areas0–2°801.50
agricultural areas/pastures211, 242, 23160sandy soil>1 m30S, W, E slopes2–5°20
Type III—agricultural/meadows and pastures areasSubtype III aagricultural areas21190loam/clay loam/silt loam>1 m45flat areas0–2°1000.60
sandy loam>1 m30
Subtype III bagricultural areas211, 242, 24375loam/clay loam/silt loam>1 m40flat areas0–2°1000.89
sandy loam>1 m30
sandy soil>1 m30
Subtype III cpastures23140peat/organic matter< 1 m30flat areas0–2°800.86
agricultural areas211, 242, 24340sandy loam>1 m20
sandy soil>1 m10
Type IV—forested areasSubtype IV a (80% forests)coniferous forest31240sandy soil>1 m50flat areas0–2°55−0.99
broad-leaved forest31120sandy loam>1 m30N, W/E slopes>5°45
mixed forest31320
Subtype IV b (90% forests)coniferous forest31240sandy loam>1 m80flat areas0–2°40−0.70
broad-leaved forest31130peat/organic matter<1 m20W/E slopesall35
mixed forest31320
Subtype IV c (90% forests)mixed forest31340sandy loam>1 m40flat areas0–2°65−0.41
broad-leaved forest31130peat/organic matter<1 m30
coniferous forest31220
Subtype IV d (95% forests)coniferous forest31285sandy loam>1 m70flat areas0–2°60−0.29
mixed forest31310sandy soil>1 m20W/E slopes2–5°15
Type V—mixed areas (with dominance of wetlands)Subtype V aconiferous forest31230peat/organic matter<1 m30flat areas0–2°80−0.15
inland marshes41125sandy soil>1 m30
pastures23120sandy soil<1 m20
Subtype V bpastures23140peat/organic matter<1 m50flat areas0–2°900.33
coniferous forest31220sandy loam>1 m15
inland marshes41120sandy soil<1 m10
sandy soil>1 m10
* % 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.
Table 12. Standardized average LST of thermal classes calculated based on Landsat 7 thermal images.
Table 12. Standardized average LST of thermal classes calculated based on Landsat 7 thermal images.
Category *Thermal Type/SubtypeMeanMinimumMaximumStandard Deviation
1I−1.91−4.36−0.300.43
2II1.50−0.214.350.49
3III a0.60−0.963.120.61
4III b0.89−1.013.130.55
5III c0.86−0.722.710.37
6IV a−0.99−2.610.520.46
7IV b−0.70−2.230.660.36
8IV c−0.40−1.711.220.39
9IV d−0.29−1.591.150.34
10V a−0.15−1.651.700.42
11V b0.33−1.312.130.37
* The first column (category) refers to 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.
Table 13. The 0/1 contingency table. When a given class from Paszyński’s classification (Figure 8) overlaps with a class from the LST classification (Figure 10), a cell takes a value of 1. Otherwise, a cell takes a value of 0.
Table 13. The 0/1 contingency table. When a given class from Paszyński’s classification (Figure 8) overlaps with a class from the LST classification (Figure 10), a cell takes a value of 1. Otherwise, a cell takes a value of 0.
Paszyński012345678910111213
LST
000000000000000
101000000000000
200111100000111
300111100000111
400111100000111
500111100000111
600000000001000
700000000010000
800000001000000
900000010100000
1000111100000111
1100111100000111
Table 14. Proportional cross tabulation: the proportion of pixels from Paszyński’s theoretical approach (Figure 8)—columns—in the LST classes; Figure 9—rows.
Table 14. Proportional cross tabulation: the proportion of pixels from Paszyński’s theoretical approach (Figure 8)—columns—in the LST classes; Figure 9—rows.
Paszyński12345678910111213
LST
10.79 *0.000.000.030.000.030.000.000.010.080.000.000.00
20.000.350.210.070.070.050.010.010.000.000.410.590.53
30.000.000.150.030.130.010.000.000.000.000.000.010.00
40.000.020.190.090.140.020.000.020.010.000.250.120.08
50.000.420.270.320.210.080.040.030.010.000.170.170.27
60.120.000.010.020.010.100.050.100.140.410.000.010.01
70.030.010.010.060.020.220.200.120.350.350.000.000.00
80.020.050.010.010.050.120.370.200.170.030.000.000.00
90.000.000.010.040.010.250.220.270.220.090.010.000.00
100.020.000.050.130.100.030.030.120.040.020.070.030.03
110.010.140.080.190.240.080.090.130.050.020.090.060.08
* Cells with high proportion values are colored: green shows the highest agreement (>0.5), followed by yellow (0.4–0.5), orange (0.3–0.4), and red (0.2–0.3). No color indicates very low or no agreement (<0.2).
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ciazela, M.; Ciazela, J. Topoclimate Mapping Using Landsat ETM+ Thermal Data: Wolin Island, Poland. Remote Sens. 2021, 13, 2712. https://doi.org/10.3390/rs13142712

AMA Style

Ciazela M, Ciazela J. Topoclimate Mapping Using Landsat ETM+ Thermal Data: Wolin Island, Poland. Remote Sensing. 2021; 13(14):2712. https://doi.org/10.3390/rs13142712

Chicago/Turabian Style

Ciazela, Marta, and Jakub Ciazela. 2021. "Topoclimate Mapping Using Landsat ETM+ Thermal Data: Wolin Island, Poland" Remote Sensing 13, no. 14: 2712. https://doi.org/10.3390/rs13142712

APA Style

Ciazela, M., & Ciazela, J. (2021). Topoclimate Mapping Using Landsat ETM+ Thermal Data: Wolin Island, Poland. Remote Sensing, 13(14), 2712. https://doi.org/10.3390/rs13142712

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop