Univerza v Ljubljani Fakulteta za gradbeništvo ingeodezijo Univerza v Ljubljani Fakulteta za gradbeništvo ingeodezijo PODIPLOMSKI ŠTUDIJ GEODEZIJE DOKTORSKI ŠTUDIJ Kandidat: KLEMEN ZAKŠEK, univ.dipl.inž.geod. DOLOČITEV TEMPERATURE ZRAKA IZ MODELA TEMPERATURE POVRŠJA VISOKE PROSTORSKE IN ČASOVNE LOČLJIVOSTI Doktorska disertacija štev.: 173 AIR TEMPERATURE DETERMINATION FROM LAND SURFACE TEMPERATURE MODEL OF HIGH SPATIAL AND TEMPORAL RESOLUTION Doktoral thesis No.: 173 Temo doktorske disertacije je odobril Senat Univerze v Ljubljani na 12. seji dne 13. februarja 2007 in imenoval mentorja izr.prof.dr. Krištofa Oštirja The theme of the doctoral thesis was approved by the Senate of the University of Ljubljana at its 12th meeting on Februar 13, 2007, who appointed Assoc. Prof. Krištof Oštir, Ph.D. as mentor. Ljubljana, 17. julij 2007 Univerza v Ljubljani Fakulteta za gradbeništvo ingeodezijo Univerza v Ljubljani Fakulteta za gradbeništvo ingeodezijo Komisijo za oceno ustreznosti teme doktorske disertacije v sestavi prof. dr. Jože Rakovec, UL FMF, dr. Tomaž Podobnikar, znanstveni sodelavec, ZRC SAZU, izr. prof. dr. Radoš Šumrada, FGG, doc. dr. Mojca Kosmatin Fras, FGG, in izr. prof. dr. Krištof Oštir, FGG in ZRC SAZU, je imenoval Senat Fakultete za gradbeništvo in geodezijo na 19. redni seji dne 20. aprila 2005. Komisijo za oceno doktorske disertacije v sestavi doc. dr. Mojca Kosmatin Fras, FGG, izr. prof. dr. Krištof Oštir, FGG in ZRC SAZU, prof. dr. Jože Rakovec, FMF UL, in dr. Thomas Holzer-Popp, Nemška vesoljska agencija (DLR-DFD-KA), je imenoval Senat Fakultete za gradbeništvo in geodezijo na 7. redni seji dne 25. aprila 2007. Komisijo za zagovor doktorske disertacije v sestavi prof. dr. Bojan Majes, dekan FGG, predsednik, doc. dr. Mojca Kosmatin Fras, FGG, izr. prof. dr. Krištof Oštir, FGG in ZRC SAZU, prof. dr. Jože Rakovec, FMF UL, in dr. Thomas Holzer-Popp, Nemška vesoljska agencija (DLR-DFD-KA), je imenoval Senat Fakultete za gradbeništvo in geodezijo na 8. redni seji dne 30. maja 2007. Univerza v Ljubljani Fakulteta za gradbeništvo ingeodezijo Univerza v Ljubljani Fakulteta za gradbeništvo ingeodezijo Commission for the adequacy of the doctoral thesis consisting of Prof. Jože Rakovec, Ph.D., UL FMF, Tomaž Podobnikar, Ph.D., scientific co-worker, ZRC SAZU, Assoc. Prof. Radoš Šumrada, Ph.D., FGG, Assist.Prof. Mojca Kosmatin Fras, Ph.D., FGG, and Assoc. Prof. Krištof Oštir, Ph.D., FGG and ZRC SAZU, was appointed by the Senate of the Faculty of Civil and Geodetic Engineering at its 19th meeting on April 20, 2005. Commission for the evaluation of the doctoral thesis consisting of Assist. Prof. Mojca Kosmatin Fras, Ph.D., Assoc. Prof. Krištof Oštir, Ph.D., Prof. Jože Rakovec, Ph.D., UL FMF, and Thomas Holzer-Popp, Ph.D., German Aerospace Centre (DLR-DFD-KA), was appointed by the Senate of the Faculty of Civil and Geodetic Engineering at its 7th meeting on April 25, 2007. Commission for the defence of the doctoral thesis consisting of Prof. Bojan Majes, Ph.D., FGG dean, president, Assist. Prof. Mojca Kosmatin Fras, Ph.D., FGG, Assoc. Prof. Krištof Oštir, Ph.D., FGG, Prof. Jože Rakovec, Ph.D., UL FMF, and Thomas Holzer-Popp, Ph.D., German Aerospace Centre (DLR-DFD-KA), was appointed by the Senate of the Faculty of Civil and Geodetic Engineering at its 8th meeting on May 30, 2007. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. IZJAVA O AVTORSTVU Podpisani KLEMEN ZAKŠEK, univ.dipl.inž.geod., izjavljam, da sem avtor doktorske disertacije z naslovom: »DOLOČITEV TEMPERATURE ZRAKA IZ MODELA TEMPERATURE POVRŠJA VISOKE PROSTORSKE IN ČASOVNE LOČLJIVOSTI«. I, the undersigned Klemen Zakšek, Univ. Geod. Civ. Eng., hereby declare that I am the author of the doctoral thesis titled: “AIR TEMPERATURE DETERMINATION FROM LAND SURFACE TEMPERATURE MODEL OF HIGH SPATIAL AND TEMPORAL RESOLUTION”. Ljubljana, 17.07.2007 ……………………………….. (podpis / signature) Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Errata Page Line Error Correction Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Bibliographic-documental information and abstract UDC 528.88:551.524/.526(043.3) Autor Klemen Zakšek, Supervisor dr. Krištof Oštir Title Air temperature determination from a land surface temperature model of high spatial and temporal resolution. Notes 143 p., 7 tab., 69 fig., 29 equ., 4 app. Key words ambient air temperature, land surface temperature, multiple regression, contextual analyses, wavelet transform, remote sensing Abstract The thesis presents the process of determining a continuous ambient air temperature (AAT) field close to the ground from the land surface temperature (LST) measured by remote sensing. Remote sensing data was used to parameterize the differences between LST and AAT by multiple regression. Most data was acquired by SEVIRI, which is an instrument on board Meteosat Second Generation. Its spatial (temporal) resolution equals approximately 5000 m (15 min) in central Europe, which is spatially too sparse to make an accurate estimation of LST or AAT in inhomogeneous areas. Therefore, SEVIRI LST was first downscaled to 1000 m spatial resolution through contextual analyses using the data about vegetation indices, albedo, elevation, incidence angle of the sun, and then to 125 m spatial resolution by wavelet transform using the simulated data of solar radiance exposure. It was discovered that the 1000 m spatial resolution provides statistically better data and AAT was then also determined in the 1000 m spatial resolution. The proposed procedure makes it possible to estimate LST and AAT in high spatial (1000 m) and temporal (15 min) resolution with the accuracy of 2.1 °C. This is the appropriate approach when AAT measurements are not available for interpolation, which usually provides results of better accuracy (1.6 °C). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Bibliografsko-dokumentacijska stran in izvleček UDK 528.88:551.524/.526(043.3) Avtor Klemen Zakšek, univ. dipl. inž. geod. Mentor doc. dr. Krištof Oštir Naslov Določitev temperature zraka iz modela temperature površja visoke prostorske in časovne ločljivosti Obseg in oprema 143 str., 7 pregl., 69 sl., 29 en., 4 pril. Ključne besede temperature zraka, temperature površja, multipla regresija, kontekstualne analize, valčna transformacija, daljinsko zaznavanje Povzetek Naloga opisuje določitev zveznega polja temperature zraka na višini 2 m iz temperature površja. Z multiplo regresijo so bile parameterizirane razlike med temperaturo površja in zraka. Večino uporabljenih podatkov je bilo zajetih s senzorjem SEVIRI na krovu satelita Meteosat Second Generation. Prostorska ločljivost teh podatkov znaša na območju centralne Evrope približno 5000 m, kar je premalo za natančno določitev temperature površja ali zraka na nehomogenih območjih. Zato je bila najprej izboljšana prostorska ločljivost temperature površja na 1000 m s kontekstualnimi analizami na osnovi podatkov o vegetacijskih indeksih, albedu, nadmorski višini in vpadnemu kotu sonca. Nadalje je bila ločljivost izboljšana še na 125 m z valčno transformacijo na osnovi simulirane osončenosti. Rezultati v ločljivost 1000 m dajejo zadovoljive rezultate, zato je bila temperatura zraka določena v isti ločljivosti. Rezultat naloge je metoda, s katero lahko dobimo temperaturo površja in zraka v visoki prostorski (1000 m) in časovni (15 min) ločljivosti s točnostjo 2,1 °C. Metoda se izkaže za ustrezno, kadar meritve temperature zraka niso na voljo za interpolacijo, ki načeloma omogoča točnejše rezultate (1,6 °C). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Acknowledgment I would like to express my gratefulness to many people and organizations that supported and helped me during my study. First of all, I would like to thank my mentor Dr. Krištof Oštir and Dr. Tomaž Podobnikar (both ZRC SAZU) for the introducing me to remote sensing and geographical information system. I also got many advices and moral support from all my colleagues at the Institute of Anthropological and Spatial Studies (ZRC SAZU). I spent ten months as a guest scientist at German aerospace centre (DLR-DFD-KA), where I worked with Marion Schrödter-Homscheidt, Dr. Thomas Holzer-Popp and Dr. Kurt Günther. I have broadened my knowledge about thermal remote sensing and meteorology under their guidance. I spent additional two months on a study visit at the Centre for urban and regional planning in Besançonu (ThéMA; France), where I studied spatial analyses under the supervision of Dr. Daniel Joly. I got the idea for the thesis by visiting Slovenian metrological office (ARSO), where I mostly cooperated with Mojca Dolinar. I was also given many useful advices by Dr. Norbert Pfeifer (TU Wien; Austria). I have to thank to all the members of the commission for the evaluation of the thesis for all their comments. I am also grateful to Sunčan Stone, Živa Zakšek and Lea Udovč for correcting errors in the manuscript. I have already mentioned many institutions that hosted me during my research; special thanks also to the Slovenian Ministry of Higher Education, Science and Technology, Ad futura foundation and French Government for financial support during my study. It was easy to research such an interesting topic because I had always had a support of my parents and friends. Therefore, thank you! Antje, you were right about me, science is indeed a game for me and I enjoy playing it. Thank you for believing in me! Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Zahvala Ob tej priložnosti bi rad izrazil svojo hvaležnost mnogim ljudem in organizacijam, katerih podporo sem bil deležen med svojim raziskovalnim delom. Najprej bi se rad zahvalil svojemu mentorju dr. Krištofu Oštirju in dr. Tomažu Podobnikarju (oba ZRC SAZU). Oba sta dala pečat mojemu študiju in me vpeljala v daljinsko zaznavanje ter geografske informacijske sisteme. Veliko nasvetov in tudi moralne podpore sem dobil od vseh sodelavcev na Inštitutu za antropološke in prostorske študije (ZRC SAZU). Deset mesecev sem kot gostujoči znanstvenik deloval v nemški vesoljski agenciji (DLRDFD- KA), kjer sem delal z Marion Schrödter-Homscheidt, dr. Thomasom Holzer-Poppom in dr. Kurtom Güntherjem. Z njihovo pomočjo sem razširil znanje na področju termalnega daljinskega zaznavanja in meteorologije. Dva meseca sem bil tudi na študijskem dopustu v centru za urbano in regionalno planiranje v Besançonu (ThéMA; Francija), kjer sem se ukvarjal s prostorskimi analizami pod mentorstvom dr. Daniela Jolya. Zamisel za obravnavano temo sem dobil na Uradu za meteorologijo (ARSO), kjer sem sodeloval predvsem z Mojco Dolinar. Veliko uporabnih nasvetov sem dobil tudi od dr. Norberta Pfeiferja z dunajske tehniške univerze (TU Wien). Za vse komentarje bi se rad zahvalil članom komisije za oceno disertacije. Hvaležen sem tudi vsem trem lektorjem: Sunčan Stone, Živa Zakšek in Lea Udovč. Omenil sem že organizacije, ki so me gostile med študijem, še posebej pa bi se rad zahvalil Ministrstvu za visoko šolstvo, znanost in tehnologijo, agenciji Ad futura in francoski vladi za finančno podporo. Raziskovalno delo mi nikoli ni pretežavno, saj sem imel vedno podporo staršev in prijateljev. Zato hvala vsem! Antje, prav si imela, znanost mi je res kot igra in uživam v njej. Hvala da mi vedno stojiš ob strani! Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Table of contents 1 INTRODUCTION 1 1.1 Structure of the thesis 2 1.2 Research background 3 1.3 Processing scheme and research hypotheses 4 2 AMBIENT AIR TEMPERATURE INTERPOLATION IN INHOMOGENEOUS AREAS 7 2.1 Methods 9 2.2 Case study: Slovenia 13 2.3 Results of the modelling procedure (Test 1) 19 2.4 Iterative improvement of the results 24 2.5 Discussion 29 3 MODELLING LAND SURFACE TEMPERATURES WITH CONTEXTUAL ANALYSES 31 3.1 Description of the LandSAF LST product 33 3.2 Correction of the SEVIRI LST 35 3.3 Downscaling 43 3.4 Case study: Slovenia 46 3.5 Discussion 57 4 MODELLING LAND SURFACE TEMPERATURES WITH WAVELET TRANSFORM 61 4.1 Discrete wavelet transform 63 4.2 Solar radiance exposure simulation 66 4.3 Case study: mountainous area on the border between Slovenia and Austria 72 4.4 Discussion 78 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 5 AMBIENT AIR TEMPERATURE PARAMETERIZATION FROM LAND SURFACE TEMPERATURE 81 5.1 Published methods 83 5.2 Testing existing methods 87 5.3 A new AAT parameterization 100 5.4 Discussion 110 6 CONCLUSION 113 7 SUMMARY 119 8 POVZETEK 123 8.1 Hipoteza 124 8.2 Interpolacija temperature zraka na nehomogenih območjih 126 8.3 Izboljšava prostorske ločljivosti z uporabo kontestualnih analiz 129 8.4 Izboljšava prostorske ločljivosti z uporabo valčne transformacije 130 8.5 Parametrizacija temperature zraka iz temperature površja 132 8.6 Zaključek 135 LITERATURE AND SOURCES 137 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. List of tables Table 1: Size of spatial analyses windows used during data preparation. 15 Table 2: Vegetation amount: reclassification (quality to quantitative) of land cover data (empirically chosen values that fit 8-bit format without units). 17 Table 3: Test 1 – frequency [%] of residuals for 7:00, 14:00 and 21:00 (365 days; 20 stations). 23 Table 4: Test 2 – frequency of residuals for 7:00, 14:00 and 21:00 (365 days; 20 stations) in comparison to Test 1. 25 Table 5: Test 3 – frequency of residuals for 7:00, 14:00 and 21:00 (365 days; 20 stations) in comparison to Test 2. 26 Table 6: : Images used to fit the SEVIRI LST to MODIS LST; there are five images with a zenith angle greater than 25°, which is approximately half of the SEVIRI zenith angle. 46 Table 7: Comparison of the characteristics between the potential methods to be used for AAT determination. 87 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. List of figures Figure 1: Processing scheme of the thesis. 6 Figure 2: Proposed procedure of determining the AAT spatial distribution. 11 Figure 3: Spatial distribution of meteorological stations (data source Environmental Agency of the Republic of Slovenia) used in the case study regarding the elevation (elevation data source: DMR 100, November 2005, © Surveying and Mapping Authority of the Republic of Slovenia). 14 Figure 4: Topographic down-up position; many exposed valleys in blue and ridges in red are present in the Alps. 16 Figure 5: Land cover data (Kokalj and Oštir, 2006). 18 Figure 6: MODIS LST data (20.05.2005 at 13:45 CET; DAAC, 2006). 19 Figure 7: Correlation index r for six explanatory variables and the seven window sizes; situation on the 01.01.2005 at 7:00. 21 Figure 8: Histogram of the residuals for 7 classes (01.01.2005 at 7:00 situation; 20 stations). 22 Figure 9: Histogram of the residuals for 7 classes (1095 situations; 20 stations). 23 Figure 10: Variation of residuals standard deviation over one year for Test 1. 24 Figure 11: Histogram of residuals for 7 classes (1095 situations; 20 stations; measurements at each hour have a unique colour). 27 Figure 12: Variation of residuals standard deviation over a year; comparison between the Test 1 and the Test 4 for measurements at 7:00 and 21:00. 28 Figure 13: Interpolated AAT field in 100 m spatial resolution; example for 04.01.2005 at 7:00. 28 Figure 14: Comparison between an independent cloud mask (APOLLO algorithm) and mask used by LandSAF for the Mediterranean area on the 02.07.2005. 34 Figure 15: SEVIRI always looks upon Europe from the south, thus it sees less shadows (in red) and a greater part of the southern area than MODIS. 36 Figure 16: The position of MSG SEVIRI in the equatorial (left) and meridian (right) plane. 38 Figure 17: The zenith angle of the satellite on the surface has a significant influence on LST. If the sensor records an image from the side (left), the image is deformed, thus the sunny areas might be over-represented compared to the shades, which is not a problem with the images recoded from the top (right). 39 Figure 18: Influence of the inclined surface upon the satellite image – a sensor detects more of the surface that is orientated towards the sensor (light coloured). 40 Figure 19: Influence of the inclined surface on the solar radiance exposure. 41 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 20: Influence of NDVI on LST [°C] during the daytime [h] (darker colours present points with a higher NDVI; 29.08.2005) on the chosen stations marked with different symbols. As expected, it can be seen that the areas covered in vegetation do not heat up as much as those with less vegetation. 44 Figure 21: Downscaling with a high spatial resolution NDVI. 45 Figure 22: Incidence angle of SEVIRI in Slovenia averaged to 1 km pixel (elevation data source: DMR 100, November 2005, © Surveying and Mapping Authority of the Republic of Slovenia). The cosine of the angle is a measurement of the surface exposure to SEVIRI. SEVIRI sees more of those areas that are more exposed to it (the incidence angle equals between 40° to 50° for most of the area). 47 Figure 23: Layers A (above) and B (below) used for land cover correction (data source © EEA, Copenhagen, 2005); look at equation (7). 49 Figure 24: Scatterplot MODIS-SEVIRI LST for the days listed in the table 6 (data for each day is in a different colour). 51 Figure 25: Scatterplot MODIS-corrected SEVIRI for the days listed in the table 6 (data for each day is in a different colour). 51 Figure 26: Difference between corrected SEVIRI and MODIS LST; an example for 07.04.2006 at 12:15 GMT. 52 Figure 27: A scatterplot MODIS-SEVIRI LST before the fit to MODIS LST; an example for 29.10.2005 at 12:15 GMT. 53 Figure 28: A scatterplot MODIS-corrected SEVIRI after the fit to MODIS LST; an example for 29.10.2005 at 12:15 GMT. 53 Figure 29: Corrected SEVIRI LST (5000 m spatial resolution); an example for 19.07.2006 at 12:15 GMT. 55 Figure 30: MODIS (at 12:20 GMT) and downscaled SEVIRI LST (at 12:15 GMT) for 19.07.2006. 56 Figure 31: A scatterplot MODIS-downscaled SEVIRI; an example for 07.04.2006. 57 Figure 32: Direct solar radiance exposure is much stronger than the diffuse radiation, thus the areas in the shade have a significantly lower LST (frost is only in the shade). 63 Figure 33: The influence of relief morphology on direct (above) and diffuse (below) solar radiance exposure. 67 Figure 34: Solar incidence angle determines whether the surface is in the hill shade. 69 Figure 35: Cast shade determination. 70 Figure 36: The sky-view factor of the cone depends on the angle between the cone and the horizontal plane. 72 Figure 37: Relief parameters that influence solar radiance exposure (Alpine area of 60 by 60 km). 74 Figure 38: The ratio between the diffuse and the direct solar radiance exposure for the 21st decade (a ten day time period) in the year; it can be seen that the ratio is greater Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. within the high mountains due to convective cloudiness typical for clear summer afternoons – the ratio of the lowland (60%) was therefore used because the sky was cloudless in the study area (Alpine area of 120 by 120 km). 75 Figure 39: LST (above) and their histograms (below) for the proxy of solar radiance exposure (left) and 1 km LST (right) within the case study area. 76 Figure 40: Downscaled (125 m) LST (left) and the reference LST (right) in the case study area. 77 Figure 41: Differences between the downscaled and reference LST. 78 Figure 42: RMSE of AAT noon forecast in Germany during the year 2002; the line with circles presents the error for 24 h AAT forecast (2.3 °C for the entire year; DWD, 2007). 82 Figure 43: TVX method: the regression line (right) parameters are determined within a moving window (light grey) for its centre cell. 84 Figure 44: Energy balance on the land surface. 86 Figure 45: DSLF for some of the chosen German meteorological stations (marked with different symbols) on 30.10.2005; notice how the values remain constant over a three hour long period. The line shows how DSLF changes during the day. 89 Figure 46: Position of the German meteorological stations used within the study – only those marked with circles and an ID number were used (data source ESRI maps and GLOBE DEM – Hastings and Dunbar, 1998). 90 Figure 47: Measured-parameterized AAT scatterplot; MODIS LST and NDVI; r = 0.93; RMSE = 2.5 °C. 92 Figure 48: Measured-parameterized AAT scatterplot; SEVIRI LST &MODIS NDVI; r = 0.96; RMSE = 2.6 °C. 92 Figure 49: Relationship NDVI-LST and EVI-LST on the 27.10.2005 south of Ljubljana (Slovenia; Zidarič, 2007). 93 Figure 50: Measured-modelled AAT scatterplot for selected days at 6:00, 13:00 and 20:00 GMT; results obtained in Slovenia with the original Meteotest parameterization; r = 0.91; RMSE = 4.0 °C. 95 Figure 51: Measured-modelled AAT scatterplot for selected days at 6:00, 13:00 and 20:00 GMT; results obtained in Slovenia with a changed Meteotest parameterization; r = 0.97; RMSE = 2.3 °C. 96 Figure 52: Measured-modelled AAT scatterplot for selected days at 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.75; RMSE = 4.5 °C; the points within the ellipse belong to the Zugspitze meteorological station. 97 Figure 53: Measured-modelled AAT scatterplot for selected days at 6:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.90; RMSE = 3.3 °C. 98 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 54: Measured-modelled AAT scatterplot for selected days at 13:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.88; RMSE = 6.4 °C. 98 Figure 55: Measured-modelled AAT scatterplot for selected days at 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.96; RMSE = 2.9 °C. 99 Figure 56: Measured-modelled AAT scatterplot for selected days at 6:00, 13:00 or 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.89; RMSE = 4.4 °C. 99 Figure 57: Measured-modelled AAT scatterplot for selected days from 4:00 to 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.90; RMSE = 3.8 °C. 100 Figure 58: Residuals (between the parameterized and measured AAT) [°C] as a function of daily time [h] (03.06.2005); each station has its own symbol. 101 Figure 59: Lines representing the residuals [°C] plotted as a function of accumulated solar irradiation [Whm-2] for selected days in 2005; more curved lines are made for winter days and less curved for summer days. 102 Figure 60: Differences between LST and AAT [°C] as a function of time of day [h] (03.06.2005); each station has its own symbol. 103 Figure 61: Lines representing the residuals [°C] plotted as a function of the time of day [h] for a few days in 2005. 103 Figure 62: Measured-modelled AAT scatterplot for selected days between 4:00 and 20:00 GMT; results obtained in Germany with a new parameterization; r = 0.96; RMSE = 2.1 °C. 104 Figure 63: RMSE variation for selected situations. 105 Figure 64: RMSE variation for selected situations (wind velocity not considered). 106 Figure 65: An AAT example for Slovenia, generated with the described procedure (27.07.2005 at 20:00 GMT) in 1000 m spatial resolution. One can observe numerous no-data areas, which are the consequence of cloud coverage. 107 Figure 66: Histogram of residuals considering wind (blue) and not considering wind (red) in the parameterization. 108 Figure 67: AAT for Slovenia (an example for 28.07.2005 at 13:00 GMT) in 1000 m spatial resolution. 109 Figure 68: Measured-modelled AAT scatterplot for selected days at 13:00 GMT; results obtained in Slovenia with the use of the improved downscaling procedure on a few stations; r = 0.97; RMSE = 1.8 °C. 110 Figure 69: Bias (above) and the accuracy estimation for both approaches of AAT determination. 114 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Abbreviations and symbols AAT ambient air temperature (measured on the height of 2 m above ground) AATSR advanced along-track scanning radiometer AL albedo ASTER advanced spaceborne thermal emission and reflection radiometer AVHRR advanced very high resolution radiometer CET Central European time DEM digital elevation model DSLF down-welling surface long-wave flux DSSF down-welling surface short-wave flux (E)TM (enhanced) thematic mapper EVI enhanced vegetation index FVC fraction of vegetation cover GIS geographic information system GMT Greenwich mean time IR infrared spectra LandSAF Land Surface Analysis Satellite Applications Facility LST land surface temperature MODIS moderate-resolution imaging spectroradiometer MSG meteosat second generation NASA National Aeronautic and Space Administration NDVI normalised differential vegetation index RMSE root mean square error SEVIRI spinning enhanced visible and infrared imager TduP topographic down-up position TIR thermal-infrared spectra TVX temperature vegetation index VIS visible spectra Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. A average angle of horizon a azimuth d horizontal distance e equation of time h height I (Moran) autocorrelation index Lc land cover morphology correction N cloud coverage (in oktas) n number of measurements v n normal vector to the surface v nE Rc normal vector to the ellipsoid relief morphology correction Rdif diffuse radiance exposure r Rdif diffuse radiance exposure influenced by the relief Rdir direct radiance exposure r Rdir direct radiance exposure influenced by the relief Rr solar radiance exposure influenced by the relief Rr proxy for the solar radiance exposure influenced by the relief proxy r (Pearson) correlation index v s vector of sunrays t time u wind speed z zenith angle W spectral radinace ß land cover roughness (average inclination of the land cover surface) . solar declination . geographical (ellipsoidal) latitude . emissivity . relief slope . incidence angle Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. . geographical (ellipsoidal) longitude µ relief aspect . ratio between diffuse and direct radiation . standard deviation . time angle . sky-view factor Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. INTRODUCTION Any description of reality is merely an abstraction of reality at a given moment – i.e. a model that simplifies reality by omitting all its insignificant aspects in order to preserve merely the important aspects for a certain purpose (Šumrada, 2006). The model structure depends on the purpose and application of the model as well as on the modelling experience. Modelling is the process of developing a model for a certain purpose. The modelling technique can be used in the analysis of complex environmental characteristics such as land surface temperature or ambient air temperature. Land surface temperature (LST) is a result of vertical energy fluxes that heat or cool the surface of the Earth. LST is becoming important in a number of research fields and applications, for example in forestry to detect forest fires, in road transport to predict frost on roads, etc. The Earth’s surface can be covered with vegetation or anthropogenic objects, thus LST, which can be detected by remote sensing techniques, is not necessarily equal to the ground temperature; however, it is correlated to it. The term ambient air temperature (AAT) presents the air temperature close to the ground. Within the thesis, only the air temperature measured on meteorological stations on the height of 2 m above ground is meant as AAT. In the broader meteorological community, one can also consider AAT as the temperature of higher atmosphere layers that can be measured by radio-sounders from the ground or satellites. AAT on the height of 2 m is important to broader community because it can influence the individual’s daily routine. Accurate AAT spatial distribution in a high spatial (1000 m) and temporal (15 min) resolution is also becoming economically interesting, because its data is needed for example by civil engineers who try to optimize buildings as regards energy use. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 1.1 Structure of the thesis The thesis describes three problems: AAT interpolation from ground based measurements, AAT parameterization from LST and LST downscaling. Each problem is fully described in its own chapter (LST downscaling is performed in two steps and described in two chapters). The introduction presents the thesis structure, the research background and the research hypotheses. Notice that “the state of the art” is not presented in the introduction but for each problem within its corresponding chapter. The interpolation of AAT measurements in inhomogeneous areas is described in the second chapter. Geostatistic and deterministic interpolation methods are based on autocorrelation, which is low in inhomogeneous areas, thus a statistical method based on multiple regression is the alternative solution in such cases. The second chapter also investigates how the spatial resolution of input data reflects on results; it is shown that the data of the higher spatial resolution generally improves the results. The described approach is optimal when AAT measurements are available at the meteorological stations. If not, using remote sensing data is an appropriate approach. Therefore, the third and fourth chapters deal with downscaling LST, which is the most important data for parameterizing AAT from remote sensing data. Two different approaches are used. The third chapter is based on contextual analyses in which the most appropriate explanatory variables are sought within a moving window and applied to the multiple regression in order to determine a higher resolution LST (5000 m input LST, 1000 m output LST). The results are compared with MODIS LST. The fourth chapter is a more theoretical one due to the lack of satisfactory data. Unfortunately, there is no appropriate ASTER LST available for the area of Slovenia, which is the optimal LST product that can be used for the validation of the downscaling procedure for spatial resolutions better than 1000 m. Therefore, a single test is made – LST is computed from the sixth channel of the Landsat TM 5 image in order to estimate the quality of the results. Wavelet transform is used to fuse the low frequency input (1000 m) LST with the high frequency solar radiance exposure data modelled in 125 m spatial resolution. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The fifth chapter discusses the transfer of LST to AAT. Three different approaches are described for two different case study areas (Slovenia and Germany). Two of these approaches were already used in the past while the third parameterization was developed a new by using multiple regression. It is shown that by using developed parameterization AAT can be determined in a calm, non-advective weather situation with a similar or even better accuracy as the meteorological forecast models. The accuracy is poorer in windy conditions. The sixth chapter presents the conclusions of the thesis. Certain specific problems are pointed out and the possibilities for future research in this area are examined. Some details (LST retrieval methods, modelling and validation methods, and case study areas description) can be found in appendixes in order to make the reading of the thesis easier. 1.2 Research background The air is not heated directly by the Sun but indirectly by the ground. Ground temperature is measured at some meteorological stations, however these measurements are not made at all meteorological stations and they are usually measured below the actual surface (depth of 5, 10, 20, 30 cm, etc.). The highest daily ground temperatures are reached in the early afternoon, while the lowest appear just before sunrise. The ground temperature (the temperature of the bulk) can differentiate from the land surface temperature (LST; the temperature of the upper millimetres of the surface skin), which can be determined with remote sensing. Today, a number of satellite sensors are capable of recording LST in medium or high spatial resolution, for example: SEVIRI (on board MSG), MODIS (on board Terra and Aqua), AVHRR (on board NOAA) and ASTER (on board Terra). These sensors operate in the visible and infrared electromagnetic spectra; images from the part of the thermal infrared spectrum can be used to generate LST. Although the spatial resolution of the listed sensors is not sufficient for some applications, satisfactory downscaling algorithms remain rare. Similar to ground temperatures, ambient air temperature (AAT) reaches its highest values in the early afternoon and the lowest just before sunrise. It is linearly correlated to elevation – it is reduced by approximately 6.5 °C for every 1000 m in the free atmosphere (Rakovec et al., Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 1998) and approximately 4 °C every 1000 m near ground in Slovenia (Mekinda-Masaron, 1995). AAT variation is spatially and temporally well predictable in homogeneous areas but the situation is worse in inhomogeneous areas (because of the relief shades, different land cover, elevation differences, etc.). It is measured at meteorological stations at 2 m above ground level, which is the standard prescribed by the World Meteorological Organization. However, the density of the meteorological network is often too sparse for high quality local temperature estimation. The installation and maintenance of new stations is usually too expensive, thus it makes greater sense to determine the AAT from other available data. GIS methods have already proven to be a convenient tool for AAT modelling: • monthly and annually average climatologic AAT fields were interpolated from AAT measurements for the entire Slovenia for the period 1961–1990 in 100 m spatial resolution (Marolt, 2000), • AAT fields for a specific situation in Slovenia were computed with a combined approach (a combination of a geostatistical and a physical model) in 1000 m spatial resolution (Dolinar, 2004), • a number of diagnostic models that make it possible to determine AAT at high spatial resolution have already been developed (Brossard et al., 2002; Joly et al., 2003; Pape and Löffler, 2004), and • if the correlation between AAT and LST is proven, then AAT can be determined from remotely sensed LST (Gallo et al., 1993; Kawashima et al., 2000; Garand et al., 2004; Meteotest, 2004; Sun et al., 2005). 1.3 Processing scheme and research hypotheses The goal of the presented thesis is to develop the methodology that can generate a continuous AAT field of a high spatial and temporal resolution: continuous AAT field can be either interpolated from measurements or parameterized from independent data. AAT was first interpolated from the stationary data in high spatial resolution. The interpolation of a continuous AAT field of high spatial and temporal resolution is possible when sufficient Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. measurements are available. In the thesis, the measurements for Slovenia (one of the case study areas) were available only for 7:00, 14:00 and 21:00 Central European Time (although most stations measure meteorological parameters every 30 min), thus only three continuous fields were possible per day. Therefore, AAT was parameterized from LST and other remote sensing data. LST can be measured with high temporal resolution only from geostationary satellites that have poor spatial resolution. Since LST is a spatially inhomogeneous parameter, it was first downscaled using contextual analyses and wavelet transform. According to the processing scheme (figure 1), four hypotheses were examined within this study. The first hypothesis is: “It is possible to model LST in a high spatial (similar or better spatial resolution as provided in the existing LST products; 1000 m or better) and temporal (15 min) resolution from geostationary satellite data (recorded by SEVIRI) by using additional data on relief and vegetation.” Second hypothesis: “It is possible to parameterize a continuous AAT field from a high spatial and temporal resolution LST model and some additional data (albedo, downwelling shortwave solar radiation flux, etc.), which can explain the difference between LST and AAT with the accuracy similar to the weather forecast accuracy that is usually between 2 and 3 °C (ARSO, 2007b; DWD, 2007).” Third hypothesis: “It is assumed that AAT is influenced merely by convection (vertical movement) and not advection (horizontal movement of air masses), because no appropriate wind field data was available for the study. Therefore, the fourth hypothesis states that the accuracy of the AAT parameterized from remote sensing data is better when the weather is calm and non-advective.” Fourth hypothesis: “When AAT measurements are available it is possible to interpolate them into a continuous high spatial resolution field by taking into consideration stationary data Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. describing the relief and vegetation with an accuracy higher than the weather forecast accuracy that is usually between 2 and 3 °C (ARSO, 2007b; DWD, 2007).” AAT measurements interpolation LST downscaling (with contextual analysis and wavelet transform) parameterization of differences between LST and AAT Figure 1: Processing scheme of the thesis. Slika 1: Diagram poteka naloge. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. AMBIENT AIR TEMPERATURE INTERPOLATION IN INHOMOGENEOUS AREAS Spatial representation of climatic data is a key point in many applications: agriculture management, environment research, transport, sport activities, etc. The climatic data has been of great importance especially during the last years due to global climate changes. An estimation of the spatial distribution of ambient air temperature (AAT) is a difficult task, because the temperature close to the ground is measured at meteorological stations (at the standard height of 2 m above ground level) that are often far apart (20–30 km). Mapping such data can be a problem since AAT can vary significantly even over small distances (even on less than 100 m). Theoretically, a measurement is strictly representative for the closest surroundings of the point where it is recorded: spatial representativity of a variable varies according to autocorrelation (the higher the autocorrelation, the broader the spatial representativity) and secondly, spatial representativity varies according to explanatory attributes. Monitored spatial temperature distribution is based on random and physical causes. The first are related to random fluctuations in measurements, specific micro scale influences, etc. The random part of information cannot be removed from the model. In contrast, physical causes, Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. which usually form spatial gradients, can be statistically described with appropriate tools. Many studies have shown that relief and land cover attributes play an important role in the spatial distribution of AAT (Anquetin et al., 1999; Benichou and Le Breton, 1987; Bolstad et al., 1998; Bootsma, 1976; Joly et al., 2003; Geiger, 1965; Yoshino, 1975). Therefore, the AAT variation as a function of elevation is a strong rule of spatial distribution. In the same way many other rules that explain AAT spatial distribution exist, thus they have to be identified and estimated. A continuous spatial field is then usually generated by interpolation from point information; various deterministic or statistic methods can be used in order to interpolate AAT. Deterministic interpolation methods are based on spatial autocorrelation: the simplest technique is inverse distance weighted (IDW) interpolation, sometimes also called the Shepard's method (Shepard, 1968). Kriging is a more complex technique based on the local spatial structure of the variable to be interpolated (Hudson and Wackernagel, 1994; Laslett, 1994; Matheron, 1963). Statistical methods introduce spatial information to the interpolation process in a different way; they are based on quantifying the correlations between the variable to be interpolated and environmental features, considered as “explanatory attributes”, such as latitude, longitude, elevation, etc. All these methods are well known and integrated within various GIS software. Furthermore, all information one needs to interpolate climatologic features can be stored and prepared within GIS (Chapman and Thornes, 2003; COST, 2000; Dobeschet al., 2001; Tveito and Schoner, 2002). The main difficulties are considering the area and the features to be interpolated and choosing a combination of methods and data to optimally describe these spatial patterns. The process is empirical: the results are evaluated step by step by removing or introducing new data to the GIS in order to find the best solution. A raster GIS model was used in this study because it easily overcomes the constraints related to spatial analysis. The presented study tends to describe the interpolation experience based on the previously stated principles. Therefore, it uses a deductive approach in order to confirm our findings and estimate the general parameters that influence the AAT spatial distribution. The defined spatial climatology requires a database of appropriate data and a set of appropriate tools. In the end, the results may be mapped or even used as a model input in some other applications. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 2.1 Methods 2.1.1. Geostatistic and deterministic interpolation methods A continuous AAT field is often generated by one of the several spatial interpolation techniques (Anderson, 2002) obtained from the AAT measured at meteorological stations. Especially kriging, a geostatistic interpolation (usually used in combination with one of the trend determination techniques), is often used for meteorological purposes. However, kriging or even deterministic interpolation (inverse distance, splines etc.), cannot be a quality interpolation technique of any variable with low spatial autocorrelation, because all these techniques are based on the assumption that objects that are close to each other are more similar than objects far apart (Kanveski and Maignan, 2004). The low autocorrelation is an important issue for geostatistic and deterministic interpolation methods, especially when the study area is inhomogeneous with a large variation of land cover and relief types between the meteorological stations. The most important positive side of the deterministic methods is their speed. On the other hand, kriging has received some criticism due to the semivariogram formulation, which requires a lot of time (Collins, 1996; Philip and Watson, 1986). Kriging also has a really strong advantage compared to the deterministic methods – it minimizes the random errors (noise) in the data to be interpolated. 2.1.2. Regression interpolation The low autocorrelation issue requires other solutions – spatial distribution of AAT can also be explained with environmental parameters/attributes of the meteorological station. These attributes can be used to explain the AAT spatial distribution if a correlation between them and the AAT is statistically significant. This means that the attributes determine the behaviour of the AAT spatial distribution. Therefore, they are first used to generate a rule, which determines the AAT behaviour as regards the explanatory attributes. Later, this rule is applied to the regularly gridded data (which contains the spatial distribution of the attributes) in order to determine the AAT spatial distribution of the entire study area. This method, which is in statistics referred to as multiple regression (see also Appendix B), can be generally used to determine the relationship between any kind of variables. The rule used to explain the Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. unknown variable with the explanatory variables, is mathematically defined as a polynomial function of explanatory attributes; the evaluated value is usually a linear combination of the explanatory attributes. Considering the geographical point of view, multiple regression seems to be a better solution for generating a continuous AAT field than the geostatistic or deterministic interpolation approach, because a deductive approach allows for true knowledge on the factors that determine AAT spatial variation. Moreover, the links between AAT and explanatory attributes highlighted by correlation reveal certain physical processes. For example, according to past experience, the link between AAT and terrain elevation, which corresponds to the physical decrease of temperature with altitude (according to the air pressure), seems to be the most significant. This method has also a drawback – a long computation time in which numerous measurements are available and a number of explanatory attributes are possible. 2.1.3. The method used The proposed statistical procedure for the determination of AAT spatial distribution, which has already been successfully used to determine a continuous AAT field (Brossard, 2002; Joly et al., 2003), consists of seven steps (figure 2): 1. computation of the AAT spatial autocorrelation (if the autocorrelation index is not close to zero, it is more sensible to use a geostatistic or deterministic interpolation instead of the proposed procedure), 2. selection of the appropriate explanatory attributes, 3. organization of the explanatory attributes data into the GIS database – projection (UTM, zone 33 N) and superimposition of the data to the common spatial resolution (100 m) and the common extent, 4. multiple regression (searching for the connection between the explanatory attributes and AAT), 5. evaluation of the results by computing the residuals and cross validation, Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 6. interpolation of the residuals (by a geostatistic or deterministic interpolation if they are autocorrelated, if not just by the trend surface), and 7. mapping the results (if the results are accurate enough, otherwise return to step 2). Selection of appropriate explanatory attributes Prepare the datasets Check the AAT autocorrelation degree (proceed if it is close to 0) Perform multiple regression for a certain situation Evaluate the results Interpolate residuals with polynomial interpolation Map the results ()yxfAATAAT ,212 += Satisfactory accuracy? YES NO ()...,, 32111 aaafAAT = 21 3 4 567 Figure 2: Proposed procedure of determining the AAT spatial distribution. Slika 2: Predlagani postopek določevanja temperature zraka pri tleh. The first three steps (figure 2) are usually preformed smoothly. The problems might arise in step number four (multiple regression), if insufficient measurements are available. Multiple regression should be performed only on those explanatory attributes that are independent and statistically significant. Therefore, the correlations between AAT and explanatory attributes are systematically estimated by the correlation index r. However, if the number of measurements is small, the computed correlation index might be a consequence of chance. Therefore, one needs to set a “line” that shows which correlation can be considered as a real correlation and not as a coincidence. Only variables that show a significant correlation index (not coincidental) are then kept for the second step, since statistically insignificant data often represents a source of errors. The significance of the correlation index was tested for all explanatory attributes using the equation following a Student t-distribution (Lowry, 2007; see also Appendix C): r · n - 2 t = 1- r 2 (1) Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. where n is the number of measurements and t is the t-distribution value. There are 18 degrees of freedom (to describe a linear relationship) when using 20 measurements, which means that if one is interested in a 95% (90%) confidence interval, the statistics t should be larger than 2.10 (1.73). In this case all correlation indexes larger than 0.38 (0.30) are statistically significant (not coincidental) and the corresponding explanatory attributes can be used within the multiple regression. The 90% confidence interval was used in the presented study. The attribute with the highest correlation index was used when there were no statistically significant attributes. All of the possible combinations of statistically significant explanatory attributes are then tested within a stepwise procedure (Efroymson, 1960). The final model corresponds to the equation that explains the most AAT variance. The fifth step (evaluation of the residuals) can also be problematic. In order to evaluate the accuracy of the results often no input data is left out in the deterministic or geostatistic interpolation method; the accuracy is then estimated by the k-fold cross or leave-one-out cross validation (Kanevski and Maignan, 2004; see also Appendix C). In the presented study the validation was especially problematic, because data for only a few meteorological stations was available (insufficient measurements to make a stable model, because if some measurements are left out of the model calibration in order to use them for the validation, the model can change significantly). Therefore, leave-one-out cross-validation was used; it involved using a single observation from the original sample as validation data, and the remaining observations as training data. This was repeated 20-times per situation so that each observation was used once as validation data – each measurement is left out once, the model is calibrated on the other 19 measurements and the standard deviation . is computed from the 20 residuals res in order to estimate the result accuracy. res .in =1 .= i (2) n -1 Furthermore, step number six is also interesting because multiple regression can already be the final step in the procedure (AAT1 is the function of attributes in figure 2). However, the results might be enhanced with additional interpolation that acts like a filter which could improve the results. In order to determine the residuals autocorrelation the residual autocorrelation index I is computed after the multiple regression is completed, because they Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. can be spatially structured at the end of the analysis even if the variable to be interpolated was not autocorrelated at the beginning of the process. If their autocorrelation level is high the residuals can be interpolated by kriging. In the opposite case (I close to 0), the trend surface is interpolated from the residuals obtained through multiple regression; the interpolation is performed as a function of the geographical position (AAT2 in figure 2) and not as a function of attributes as is the case in multiple regression (AAT1 in figure 2). On the end, when both, multiple regression and the additional regression are preformed, the results can finally be mapped. 2.2 Case study: Slovenia The proposed procedure was tested in Slovenia which covers approximately 20 000 km2. Its topography features a small coastal strip on the Adriatic Sea, the Alpine region in the northwest, hills in the central part and the Panonian plain in the east. Due to its position on the crossing of various climatic regions, Slovenia has a number of different climate types. In order to describe this variety one needs various data: • Environmental Agency of the Republic of Slovenia provided the AAT measurements, • Surveying and Mapping Authority of the Republic of Slovenia provided the digital elevation model, • land cover data was prepared by Kokalj and Oštir (2006), and • MODIS satellite data was acquired from NASA. 2.2.1. Climatic temperature data AAT measurements (at the height of 2 m) for 20 climatic stations in Slovenia were at disposal for the presented study. The chosen case study dataset covers the entire year 2005 with AAT measurements three times a day: 7:00, 14:00 and 21:00 Central European Time (CET); this means that 1095 situations (365 days times 3 daily measurements) had to be interpolated. In Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. order to have an idea as regards the interpolation quality during the year, the period was cut into four seasons: winter, spring, summer and autumn. The network of the available stations does not have an optimal geometrical distribution, because some of the stations are quite close to each other and have no significant differences regarding their attributes. The spatial distribution is also not optimal as regards elevation, for there are only five stations over 500 m and not a single station over 1000 m above sea level. The spatial and elevation distribution can be seen on the elevation map (figure 3); approximately 10% of Slovenia is at least 1000 m above sea level, but unfortunately there were no measurements from such an elevation available for our study. Figure 3: Spatial distribution of meteorological stations (data source Environmental Agency of the Republic of Slovenia) used in the case study regarding the elevation (elevation data source: DMR 100, November 2005, © Surveying and Mapping Authority of the Republic of Slovenia). Slika 3: Prostorska razporeditev uporabljenih meteoroloških postaj (vir podatkov Agencija Republike Slovenije za okolje) glede na višino (vir podatkov DMR 100, november 2005, © Geodetska uprava Republike Slovenije). The range of elevation or any other explanatory attributes on the meteorological stations is important for the quality of the results – for example, one cannot expect high accuracy results Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. for an area in the high mountains, because AAT measurements were not available for the area with such attributes. Meteorological stations are usually situated on a relatively flat relief, thus the steepest slope of an area with a meteorological station has no greater inclination than 16° (at the spatial resolution of 100 m) and this represents another problem. Slopes with an inclination of 30° or more are frequent in Slovenia, thus AAT for this case had to be determined by extrapolation, and this always brings a high level of uncertainty – the accuracy of extrapolation decreases as the dissimilarity between the attributes on the meteorological stations and the attributes in an arbitrary point rises. Therefore, it was decided that AAT spatial distribution is computed only for areas with attributes similar to the attributes of the meteorological stations; an area is “valid” when none of its significant attributes are a certain value (10% of the attributes range) below the minimum or above the maximum of the attribute’s value. 2.2.2. Derivation of explanatory variables Explanatory variables were derived from the digital elevation model (DEM), land cover data and the MODIS satellite data: Derivates of DEM Derivates of land cover MODIS satellite data elevation artificial vegetation index normalized differential vegetation index relief slope distance from the sea enhanced vegetation index relief aspect distance to the forest land surface temperature quasi-global solar irradiation relief roughness topographical down-up position Some of this data has a micro-local and some has a regional influence. Therefore, some explanatory variables used in this study were estimated within spatial analyses windows of seven different sizes centred on each meteorological station (table 1). The various sizes of spatial analyses windows made it possible to consider the effect of the surrounding area on the measured AAT at different scale levels. Table 1: Size of spatial analyses windows used during data preparation. Preglednica 1: Velikost oken analize uporabljenih med pripravo podatkov. Window number 1 2 3 4 5 6 7 Size [pixel] 3 11 29 51 101 151 201 Size [km] 0.3 1.1 2.9 5.1 10.1 15.1 20.1 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 2.2.2.1 Digital elevation model As the depended variable, AAT had to be explained with explanatory variables that were presented by different datasets describing relief and land cover characteristics. The source for the relief was a DEM with a spatial resolution of 100 m (Podobnikar and Mlinar, 2006; figure 3). Elevation, relief slope (elevation gradient), relief aspect (direction of relief slope), quasi- global solar irradiation (solar radiance exposure received by the surface on a typical cloudless day, computed to the equinox) day and topographic roughness (Gardner et al., 1990), were estimated from DEM within the seven different size windows centred on each meteorological station (table 1). Figure 4: Topographic down-up position; many exposed valleys in blue and ridges in red are present in the Alps. Slika 4: Sloj, ki opisuje višinsko razmerje med posamezno točko in okolico; na območju Alp je vidnih mnogo izpostavljenih dolin v modri in grebenov v rdeči barvi. A sixth variable, called topographic down-up position (TduP; figure 4), which is the mean elevation difference between a point and its surrounding area, was also derived from the DEM, however, it was estimated only for the 11 × 11 pixels window. A TduP negative value Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. means that the chosen point is surrounded by higher areas and the opposite holds true (surrounding areas of the chosen point are lower than the point) if TduP is positive. The spatial pattern of TduP is similar to the pattern of the relief curvature, which shows concave and convex areas. AAT is generally well correlated to TduP (with a positive correlation index r): cold air flows down to the deepest places while the crests (dominant places) are well ventilated. 2.2.2.2 Land cover data The original land cover data (Kokalj and Oštir, 2006; figure 5) generated from Landsat ETM+ images has a 25 m spatial resolution. Land cover data is especially important because it contains the information on vegetation status, which influences the evapotranspiration. The vegetation status can also be described by the normalized difference vegetation index (NDVI); the most suitable NDVI data would be a year composite (if one works only with stationary data) of a spatial resolution of 100 m or below, however no appropriate data was available. Because it was assumed that the spatial variety of the land cover in Slovenia has a large effect on AAT, a new vegetation index was generated by assigning new values to the land cover classified image regarding the average NDVI for every class (table 2) – a transformation from qualitative (land cover classes) to quantitative (NDVI) data. In order to consider the data at different scale levels, the vegetation quantity was later on also averaged like DEM derivates according to the seven analysis window sizes. Table 2: Vegetation amount: reclassification (quality to quantitative) of land cover data (empirically chosen values that fit 8-bit format without units). Preglednica 2: Količina rastja: reklasifikacija vrste rabe tal iz kvalitatinega v kvanitiativni sloj (empirično določene vrednosti, ki ustrezajo 8-bitnemu zapisu). Class number Land cover type Vegetation amount 1 urban 5 2 high building density 50 3 low building density 100 4 coniferous forest 200 5 deciduous forest 250 6 mixed forest 225 7 bush 150 8 water 125 9 meadow and filed 125 10 rock 5 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 5: Land cover data (Kokalj and Oštir, 2006). Slika 5: Vrsta rabe tal (Kokalj in Oštir, 2006). Two additional parameters were derived from land cover data – distance to the forest, which is an important parameter because air usually seems cooler within or next to a forest as in the other areas because of evapotranspiration. The second parameter derived from land cover data was the distance from the sea, which is important because of the influence of the Mediterranean climate, which is generally different from the inland climate because of the sea (AAT is lower in spring, greater in autumn). The “window process” was not applied for the last two variables. 2.2.2.3 MODIS satellite data NASA provides various data that influences AAT on the web (DAAC, 2006). The Moderate Resolution Imaging Spectroradiometer (MODIS) on board Terra and Aqua satellites provides 16-day composite vegetation indices in three different spatial resolutions containing NDVI and EVI – enhanced vegetation index. Both indices at 1000 m spatial resolution (MOD13A2 product; DAAC, 2006) were tested for three random days within the year 2005 (25.05.2005, Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 12.07.2005 and 15.10.2005). Both indices at 250 m spatial resolution (MOD13Q1 product; DAAC, 2006) were also tested on 25.05.2005. MODIS is also an excellent source of land surface temperature (LST), which is highly correlated to AAT. The LST in 1000 m spatial resolution data (MOD11L2 product; DAAC, 2006) for 20.05.2005 at 13:45 CET was used as a sample of the LST spatial distribution (figure 6). This data changes over time and is not stationary data as, for example, elevation or slope, but it can provide information, as to which areas are in general warm and which cold. Figure 6: MODIS LST data (20.05.2005 at 13:45 CET; DAAC, 2006). Slika 6: MODIS-ova temperatura površja (20. 05. 2005 ob 13:45 po srednjeevropske času; DAAC, 2006). 2.3 Results of the modelling procedure (Test 1) The modelling was applied for each of the 1095 available temperature situations. First, the autocorrelation of AAT was determined using the autocorrelation index I, which was close to zero in all the situations. This means that the spatial distribution of AAT in Slovenia is random from the statistical point of view (when only spatial distances among meteorological Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. stations are considered). Therefore, in order to determine AAT spatial distribution, multiple regression was confirmed to be more suitable than any interpolation method. The used procedure can be considered as an iterative process, in which explanatory variables can be added or removed from the analysis in order to improve the results. The influence on the interpolation of different explanatory variables has been tested before satisfactory final results were achieved. Not all described explanatory variables were prepared already during Test 1 – MODIS satellite data was not used and the distance from the sea was not considered at the time. Test 1 consisted of computing the linear correlation between AAT data and the six explanatory variables derived from DEM (elevation, slope, aspect, theoretical quasi-global radiation energy, topographic roughness and TduP) as well as the two explanatory variables derived from land cover data (vegetation quantity and distance to forest). Each of these eight variables (except TduP and the distance to the forest) was processed within the frame of seven adaptive windows in order to discover (further in the modelling procedure) the most significant spatial pattern. In the end 44 different layers were considered in the modelling procedure as factors able to explain the temperature distribution. 2.3.1. Correlation indexes for the 1st of January at 7:00 situation The example for the first set of measurements in the AAT dataset (01.01.2005 at 7:00) was used to show how correlation indexes between AAT and explanatory variables vary with scale. Figure 7 shows the correlation index r values for different window sizes. The correlation index values for the two additional factors (distance to the forest and TudP) are respectively 0.33 and 0.45. The elevation effect remains constant from the first to the seventh analysis window (the curve is almost flat), reflecting that the correlation between the elevation and AAT changes gradually in the case study area (non-scalar pattern). The correlation between AAT and elevation measures the intensity of the vertical thermal gradient. These calculations show that the elevation induces a uniform influence on this gradient (the higher the elevation, the lower the AAT); the correlations remain identical and stable across the range of the seven windows. On the other hand, the remaining explanatory attributes exhibit a strong scale effect. For topographical roughness and quasi-global radiation energy, the optimal results (the strongest Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. correlation index) are obtained with windows 4 (5 by 5 km; table 1); for vegetation quantity with window 5 (10 by 10 km); for aspect and for gradient with size 1 (0.3 by 0.3 km). The scale also affects solar radiance exposure and topographical roughness, which are the least correlated with AAT in this situation. The scale of the explanatory variables indicates the hierarchy of variables that can explain the AAT distribution in a random situation. The optimal window size varies from one situation to another and might differ completely from those in the example presented for 01.01.2005. When the conditions are similar for the entire study area (for example on a clear night), large window sizes are more suitable for the analysis. Figure 7: Correlation index r for six explanatory variables and the seven window sizes; situation on the 01.01.2005 at 7:00. Slika 7: Korelacijski koeficient za šest pojasnjevalnih spremenljivk in sedem velikosti oken analize; situacija za 01. 01. 2005 ob 7:00. Among the eight variables introduced in the systematic linear correlation, only five of them were characterised by a significant correlation index and then used within the multiple regression: • elevation (window 5; r = -0.58), • relief roughness (window 1; r = 0.47), • vegetation quantity (window 5; r = -0.47), Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. • TudP (window 2; r = 0.45), and • relief slope (window 1; r = -0.42). The listed explanatory variables were used in order to determine the AAT for the chosen situation. The standard deviation of residuals for the 01.01.2005 at 7:00 situation is 2.2 °C. The residuals can be quite high: the minimum is -4 °C (station 102) and the maximum is 4.5 °C (station 51). As expected, the mean value belongs to the class with residuals between -1 °C and 1 °C (figure 8). 60 50 40 30 20 10 0 below-5 -5to -2.5 -2.5to -1 -1 to 1 1to 2.5 2.5to 5 above 5 Residuals [°C] Figure 8: Histogram of the residuals for 7 classes (01.01.2005 at 7:00 situation; 20 stations). Slika 8: Histogram odstopanj (01. 01. 2005 ob 7:00; 20 postaj). 2.3.2. Residuals for the 1095 situations of the whole year 2005 All of the remaining 1094 situations for the year 2005 are systematically analysed in the same way as the previous example (table 3). The standard deviation of the residuals for the entire 1095 situations equals 1.8 °C. The minimum and the maximum residuals are respectively -14 °C (station 144, 06.07.2005, 7:00) and 15.5 °C (station 144, 02.03.2005, 7:00). Twelve stations have at least one (absolute) residual higher than 5 °C; among these, some have residuals lower than -7 °C (stations 97, 38, 102, 359) or higher than 7 °C (stations 464, 51, 192). Over 50% of the residuals can be found between -1 °C and 1 °C (figure 9). The 7:00 situations have higher residuals than the other situations (table 3). Frequency [%] 35 20 20 15 10 0 0 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 0.9 5.8 16.7 50.7 18.8 6.2 0.9 0 10 20 30 40 50 60 below -5 -5 to -2.5 -2.5 to -1 -1 to 1 1 to 2.5 2.5 to 5 abov e 5 Residuals [°C] Frequency [%] Figure 9: Histogram of the residuals for 7 classes (1095 situations; 20 stations). Slika 9: Histogram odstopanj (1095 meritev; 20 postaj). Table 3: Test 1 – frequency [%] of residuals for 7:00, 14:00 and 21:00 (365 days; 20 stations). Preglednica 3: Test 1 – frekvenca odstopanj ob 7:00, 14:00 (365 dni; 20 postaj). Temperature class [°C] 7:00 14:00 21:00 <-10 0.1 0.0 0.0 -10 to -5 1.4 0.5 0.6 -5 to -2.5 7.0 5.6 4.8 -2.5 to -1 15.9 17.8 16.5 -1 to 1 48.3 50.0 53.5 1 to 2.5 19.2 18.6 18.6 2.5 to 5 6.9 6.3 5.4 > 5 to 10 1.1 1.1 0.5 >10 0.1 0.0 0.0 In the same way, Test 1 shows some large residuals; especially during the first three months of the year (marked as winter in figure 10). The autocorrelation index I of the residuals was close to zero, thus a trend surface interpolated from the residuals obtained by the multiple regression was added to the residuals in order to improve the statistics. However, the results were still not satisfactory due to some extreme residuals. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 10: Variation of residuals standard deviation over one year for Test 1. Slika 10: Variabilnost standardnega odklona odstopanj prek celega leta za Test 1. 2.4 Iterative improvement of the results One could consider Test 1 as unsatisfactory due to the high number of high residuals, thus three additional solutions were tested. The procedure started again at step two – the selection of appropriate explanatory variables (figure 2). 2.4.1. First improvement: introducing distance to the sea (Test 2) The largest residuals of the Test 1 affected the stations 97, 102 and 464, which are located close to the sea. Therefore, an additional attribute was introduced into the analysis – distance to the sea. All the 44 layers introduced in the Test 1 as explanatory variable are kept; with the new one, all the 45 variables were correlated with AAT. The additional data improved the accuracy on stations 97 and 102 that were 10 to 20 km away from the sea, but the high residuals remained in the stations 464 just on the seacoast. Additional meteorological measurements were checked in order to explain the residuals but the recorded meteorological measurements on the problematic station did not significantly differentiate from the other stations. The problem was the land cover data; the meteorological station is situated just next to a saltpan, thus the area was classified as the sea (50 m from the coast), which caused a systematic error in the multiple regression. The problem was solved by “moving” the station Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 200 m to the east in order to “place” it on solid ground; afterwards smaller residuals were noticed. Test 2 results are more accurate then Test 1 results (table 4); the global standard deviation has slightly improved to 1.7 °C (1.8 °C in Test 1). At 14:00 the situation improves significantly. The minimum and maximum residuals are respectively -13.6 °C station 359, 10.02.2005 at 7:00) and 12.9 °C (station 18, 10.02.2005 at 7:00). The residuals are still extremely low, yet they are slightly improved for stations close to the sea: station 102 (-9.7 °C vs. -12.7 °C in Test 1) or 464 (4.8 °C vs. 7.5 °C in Test 1) and station 321 which is far from the sea (2.6 °C vs. 5.2 °C in Test 1). Table 4: Test 2 – frequency of residuals for 7:00, 14:00 and 21:00 (365 days; 20 stations) in comparison to Test 1. Preglednica 4: Test 2 – frekvenca odstopanj ob 7:00, 14:00 (365 dni; 20 postaj v primerjavi s Testom 1. Temperature class Frequency [%] at 7:00 Frequency [%] at 14:00 Frequency [%] at 21:00 [°C] Test 2 Test 1 Test 2 Test 1 Test 2 Test 1 <-10 0.1 0.1 0.0 0.0 0.0 0.0 -10 to -5 1.1 1.4 0.3 0.5 0.5 0.6 5 to 2.5 6.7 7.0 4.8 5.6 5.1 4.8 2.5 to 1 15.9 15.9 18.0 17.8 16.1 16.5 1 to 1 49.8 48.3 52.2 50.0 53.2 53.5 1 to 2.5 19.0 19.2 18.8 18.6 19.2 18.6 2.5 to 5 6.7 6.9 5.1 6.3 5.4 5.4 5 to 10 0.8 1.1 0.8 1.1 0.4 0.5 >10 0.1 0.1 0.0 0.0 0.0 0.0 2.4.2. Second improvement: removing distance from forest (Test 3) The problem with large residuals at station 464 appeared only in situations in which the distance from the forest was one of the explanatory attributes. Only coniferous, deciduous and mixed forests were at first considered as forests while areas covered with bush were not taken into consideration. However, areas classified as bush (mostly situated in the costal and Karst area) are actually covered with a combination of low pine forest and bush, thus the bush was reclassified as forest. The results improved with the change in the distance from the forest data. Test 3 considered removing the distance from the forest from the analysis. The results improved in comparison to the case in which the bush was not used and they were insignificantly worse when compared to the case in which the bush was used as a part of the forest. Therefore, it was decided that the distance from the forest as introduced in the modelling represents an insignificant parameter in the spatial distribution of AAT in Slovenia. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Test 3 (removed distance from forest) once more improved the results (table 5); the global standard deviation is 1.66 °C (1.82 °C, 1.72 °C in Test 1 and Test 2). The improvements are the greatest for the situations at 14:00 and 21:00. The minimum and maximum residuals are -10.5 °C (station 359, 09.02.2005 at 7:00) and 11.0 °C (station 174, 10.02.2005 at 7:00). In comparison to Test 2, the residuals are lower for stations 464 (-1.3 °C vs. 4.8 °C in Test 2), 102 (-8.3 °C vs. -9.7 °C in Test 1) and 38 (-4.5 °C vs. -5.5 °C in Test 1). Table 5: Test 3 – frequency of residuals for 7:00, 14:00 and 21:00 (365 days; 20 stations) in comparison to Test 2. Preglednica 5: Test 3 – frekvenca odstopanj ob 7:00, 14:00 (365 dni; 20 postaj) v primerjavi s Testom 2. Temperature class Frequency [%] at 7:00 Frequency [%] at 14:00 Frequency [%] at 21:00 [°C] Test 3 Test 2 Test 3 Test 2 Test 3 Test 2 <-10 0.0 0.1 0.0 0.0 0.0 0.0 -10 to -5 1.0 1.1 0.4 0.3 0.4 0.5 5 to 2.5 6.3 6.7 4.2 4.8 4.6 5.1 2.5 to 1 16.4 15.9 17.0 18.0 17.0 16.1 1 to 1 49.9 49.8 54.0 52.2 54.0 53.2 1 to 2.5 19.0 19.0 19.3 18.8 19.3 19.2 2.5 to 5 6.5 6.7 4.4 5.1 4.4 5.4 > 5 to 10 0.8 0.8 0.4 0.8 0.3 0.4 >10 0.0 0.1 0.0 0.0 0.0 0.0 2.4.3. Third improvement: introducing MODIS data (Test 4) The MODIS products were introduced into the analysis in Test 4 in order to check the accuracy of the artificially produced vegetation index from the land cover data and to estimate the correlation between LST with AAT. It was proved that LST is usually a significant attribute for AAT distribution since the r was usually higher than 0.5; it should also be taken into consideration that the effect of LST varies through the day and the year, thus the correlation of the appropriate data (for the appropriate date and time) should be better correlated. The effect of MODIS NDVI (1000 m spatial resolution) was compared to the effect of the artificially produced vegetation index on AAT. As expected, the correlation coefficients for the MODIS NDVI were closest to the correlation coefficients of the artificially produced vegetation index in window 2 due to their similar spatial resolution (table 1). Test 4 was performed once again with MODIS NDVI with 250 m spatial resolution because in most cases the gain of a spatial resolution improves the interpolation. The results were in fact Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. slightly better than when using NDVI with 1000 m spatial resolution. In the final model (used LST with 1000 m and NDVI with 250 m spatial resolution), the highest residuals remained stable but the frequency of residuals close to 0 slightly increased and the standard deviation also improved to 1.60 °C (1.66 °C after Test 3). Figure 11 presents the histogram of residuals of the final results (Test 4) according to 3 sets of measurements (for each hour). Frequency [%] 0.0 10.0 20.0 30.0 40.0 50.0 below -5 -5 to-2.5 above 5 14:00 21:00 -2.5 to-1 -1 to 1 1 to 2.5 2.5 to 5 Residuals [°C] 7:00 Figure 11: Histogram of residuals for 7 classes (1095 situations; 20 stations; measurements at each hour have a unique colour). Slika 11: Histogram odstopanj (1095 meritev; 20 postaj; meritve za posamezno uro so obarvane s svojo barvo). 2.4.4. Global quality improving and final results Figure 12 shows how the quality of the results improved through the described procedure of selecting appropriate explanatory variables. At 7:00, standard deviation improved by 0.2 °C (winter and autumn) or 0.3 °C (spring and summer). At 21:00, (and the same at 14:00) the improvement was slightly lower (0.1 °C to 0.2 °C). The accuracy remained lowest during the winter and highest during the summer. Therefore, additional improvement might be possible after reconsidering some additional explanatory variables that might explain the AAT spatial distribution during the winter. A continuous AAT field, which is the final result, was obtained in the final step. One can see (an example in the figure 13) that there are some no-data areas, because of the problems with the extrapolation (described in the description of AAT measurements). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 12: Variation of residuals standard deviation over a year; comparison between the Test 1 and the Test 4 for measurements at 7:00 and 21:00. Slika 12: Variabilnost standardnega odklona odstopanj prek celega leta; primerjava med Testom 1 in 4 za meritve ob 7:00 in 21:00. Figure 13: Interpolated AAT field in 100 m spatial resolution; example for 04.01.2005 at 7:00. Slika 13: Interpolirana temperatura zraka v prostorski ločljivosti 100 m; primer za 04. 01. 2005 ob 7:00. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 2.5 Discussion The results of the presented study, which were not validated with independent data because of the low number of measurements (leave one out cross validation was therefore used), can be compared with previous studies. Anderson (2002) tested different interpolation methods in the Phoenix metropolitan area, where the density of measurements was approximately ten-times larger than in the presented study; the best results were obtained by kriging interpolation (RMSE = 1.62 °C). The presented results are also better than the 24 hour weather forecast, but they cannot be compared to the results of the AAT high resolution interpolation on Svalbard (Joly et al., 2003), where AAT was modelled with the accuracy of 0.3 °C based on 50 measurement stations covering the area of only 8 km2. One of the goals of this study was to demonstrate the possibilities for multiple regression. It is widely used in the GIS community because it can be extremely helpful in explaining processes that are analytically difficult to describe. The results of the multiple regression analysis can provide high quality results when there are enough measurements and explanatory data at disposal. The presented case study has shown that data from a greater number of meteorological stations is necessary in order to obtain a high accuracy AAT field – more measurements would increase the reliability of the results and this would make it possible to introduce the polynomial of the second (quadratic form) or third level (cubic form) into the multiple regression (linear multiple regression was used in the study). Furthermore, more measurements mean smaller distances between meteorological stations and if their density is high enough, the autocorrelation within measurements might rise to a level where another interpolation technique might be used as an alternative method; in this case the best one would probably be co-kriginig, which also considers the correlation with other explanatory variables as weights in the interpolation. The second study goal was to find the most significant attributes that can explain the AAT spatial distribution. As expected, elevation, NDVI and LST are often among the significant attributes. Either solar radiance exposure or relief aspect are also usually significant during daytime in clear-sky conditions. NDVI and elevation are the most significant in cloudy conditions. Therefore, many attributes are correlated to AAT but the connection between Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. these attributes and the measured AAT is different for different weather situations. Furthermore, the scale (the size of the analysis window used to determine the attribute) also has an effect on the results, which makes AAT interpolation a difficult task. The described procedure uses only stationary (DEM and land cover data with their derivates prepared within different window sizes) and some data that was treated as representative for the entire year although it corresponds only to a specific situation (for example quasi-global solar radiance exposure on an equinox). The procedure offers a chance to simply add or remove the explanatory variables if one is not satisfied with the results (four different tests were done within this study). The introduction of the new explanatory variable is easy if no new measurements are required. It was shown in this chapter that a continuous AAT field can be successfully interpolated from AAT measurements. However, these measurements are not always available (perhaps there are no measurements at all in some regions – they might be too expensive, too seldom, etc.). It was decided to parameterize AAT from remote sensing data, because it was shown in chapter 2.4.3 that AAT is well correlated with LST. LST is measured in a reasonably high resolution (1 km) by polar orbiting satellites that usually make two overpasses per day. An alternative solution could be the use of geostationary satellites, which have much better temporal resolution (an hour or even less – typically every 15 minutes) but worse spatial resolution. Due to the fast weather changes in the inhomogeneous areas it would be most desirable to work with data of high temporal and spatial resolution. This means that a new GIS has to be built, one that could also use dynamic data. The first objective of the new system is to downscale the LST obtained by the geostationary satellite to the spatial resolution comparable to the LST resolution obtained by the polar orbiting satellites (described in the next chapter). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. MODELLING LAND SURFACE TEMPERATURES WITH CONTEXTUAL ANALYSES Land surface temperature (LST) is an important parameter in a variety of climatic, hydrologic, agricultural, biogeochemical, and change detection studies. It depends on net radiation, the sensible heat flux, the latent heat flux and the ground heat flux. LST is an excellent indicator of the energy balance on the Earth's surface, thus it became (on the global scale) an indicator of global warming (Dash, 2004). On the other hand, an average LST composite of a low spatial resolution is inappropriate for other applications – for example, high spatial and temporal resolution data can be important in agricultural applications because surface temperature can change fast; it was shown already by Ansari and Loomis (1959) that a pepper leaf can get warmer by 10°C in less than two minutes. LST is also a spatially very heterogeneous natural variable; shady areas can be much colder than the areas exposed to the sun. Moreover, even areas that have been exposed to the Sun for the same period of time can have a different LST due to the material characteristics such as heat capacity (the measure of energy increase per degree of temperature rise), thermal conductivity (the rate at which the heat passes through a specific thickness of a material), thermal inertia (resistance to Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. temperature change) and emissivity (ratio of energy radiated by the material to energy radiated by a black body). A number of satellite sensors operate in the thermal infrared spectrum but none is truly capable of retrieving LST in a high temporal and spatial resolution. For example, one can use Landsat ETM+ images in order to derive LST with a high spatial resolution (60 m) but its temporal resolution is extremely poor (usually only a few images of the same area per year). One can of course choose a better temporal resolution such as provided by the MSG SEVIRI images but in this case the spatial resolution is relatively low – 3 km in nadir. The compromise is usually done by using polar orbiting satellites that have a revisiting time of approximately one day; this means that a day and a night image can be retrieved every day (usually the satellite passes the area in the ascending node during the day and in the descending node during the night), but even this might be too rare for some precise Soil- Vegetation-Atmosphere transfer models. More about LST retrieval is written in Appendix A. An alternative option is modelling – one can combine high spatial resolution LST data with the diurnal cycle derived from high temporal resolution LST data or one can use high temporal resolution LST data and downscale it by integrating the spatial distribution of LST correlated data. Both approaches have already been tested. Sun and Pinker (2005) implemented the GOES-based LST diurnal cycle (30 min temporal resolution) to AVHRR (1100 m spatial resolution). LST determined from real-time GOES-8 observations has a RMSE of about 2°C and the proposed fitting algorithm provides LST with RMSE similar to the LST retrieved directly by GOES-8 but in a higher spatial resolution. A similar approach was also tested by Jin and Dickinson already in 1999, but they used GOES-8 data merely for validation. Venkateshwarlu et al. (2004) tested the other approach by applying neural networks to compute LST from old Landsat TM (not ETM+) images. They also used the data provided from non-thermal channels (LST is otherwise computed from the TIR channel of 120 m spatial resolution) because of their higher spatial resolution (30 m). The presented study is based on LST data retrieved by the Spinning Enhanced Visible and Infra-Red Imager (SEVIRI) on board the Meteosat Second Generation (MSG) because of its excellent temporal resolution; they are processed by the Land Surface Analysis Satellite Applications Facility (LandSAF, 2007). In central Europe – where the case study (Slovenia) Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. is situated – SEVIRI LST has approximately 5000 m spatial resolution. The data spatial resolution is relatively low but its temporal resolution of 15 min is excellent. The comparison with MODIS LST shows that the SEVIRI LST is greater, thus a fit of SEVIRI to the MODIS LST was done first, since MODIS LST product already contains a large archive with validated data. After the correction of the 5000 m SEVIRI LST, additional data (NDVI, digital elevation model and land cover) was used in a downscaling procedure to 1000 m based on a contextual spatial analysis. The downscaling is a difficult task because of lack of adequate data sources for the results validation – only MODIS and AVHRR sensors are appropriate for the validation, because they are not on board of only one satellite (MODIS on Terra and Aqua, AVHRR on NOAA). If a sensor is on board a single satellite, it usually crosses the area once per day, which is not enough to describe the LST daily cycle. If the sensor is on two or more satellites and it also returns during night time, more data can be used for the validation. 3.1 Description of the LandSAF LST product The LandSAF (2007) LST product is generated from the data retrieved by the SEVIRI. An important parameter for LST retrieval by remote sensing is land surface emissivity (.), which is independently estimated as a function of a fraction of vegetation cover (FVC) and land cover classification from the used data. The LST retrieval is based on the measurements from the SEVIRI channel 9 (IR 10.8 µm) and 10 (IR 12.0 µm). LST is measured every 15 min, thus the values can be determined 96 times per day. LST is derived from the measurements by the generalised split-window algorithm (Wan and Dozier, 1996) in cloud-free pixels only. The identification of pixels contaminated by a cloud, the temperature of which is much lower than the surface temperature (hence the cloud mask is a significant parameter for the LST quality), is based on the cloud mask generated by the Nowcasting and Very Short Range Forecasting Satellite Application Facility (SAFNWC) software. A comparison between a LandSAF derived cloud mask and an independently derived mask from the SEVIRI data obtained through the use of the APOLLO algorithm (Kriebel et al., 2003) showed certain differences, however the cloud mask used by LandSAF LST usually removes more clouds then necessary (figure 14), thus one can rely there are no gross errors due to the bad cloud mask (even though Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. some sub-pixel clouds might be found – for example on MODIS images). The quality of the derived LST also depends on the sensor performance, the accuracy of atmospheric corrections and the spectral variation in emissions from different land-surface elements. Another significant parameter of the LST quality is land surface heterogeneity that can produce a significant variation in LST measurements as a function of the view zenith angle. The expected LST accuracy is approximately 2 °C (the possibility of large errors for large zenith angles or wet conditions). Figure 14: Comparison between an independent cloud mask (APOLLO algorithm) and mask used by LandSAF for the Mediterranean area on the 02.07.2005. Slika 14: Primerjava med neodvisno masko oblakov (APOLLO) in masko uporabljeno pri LandSAFu na območju Sredozemlja za 02. 07. 2005. The LandSAF LST product is still in its pre-operational phase, which means that it has already been scientifically checked but not validated. So far SEVIRI LST has been compared with MODIS and AATSR LST. Madeira et al. (2005) compared SEVIRI LST with LST from both mentioned sensors. They suggested that the differences between various LST data might occur due to the different emissivity data used by MODIS and AATSR, undetected clouds, or cloud shadows affecting only one of the products, different zenith angles, and thus different perspectives of similar scenes. They have discovered that MODIS LST is always 2–5°C lower while AATSR LST is always a bit greater than SEVIRI LST. Noyes et al. (2006) made their comparison only with AATSR LST. They confirmed the results from the previous study by estimating that AATSR LST is generally 1–2°C greater than SEVIRI LST. Their results also Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. indicated that the observed biases between the AATSR and SEVIRI LST data are seasonal, wand that the maximal difference (SEVIRI minus AATSR) occurs during the summer months. 3.2 Correction of the SEVIRI LST The comparisons between SEVIRI, MODIS and AATSR LST data showed some systematic differences that cannot be neglected in the downscaling procedure. Therefore, the mean difference has to be minimized and the seasonal bias has to be removed before the downscaling procedure starts. The MODIS LST is a standard MODIS product (MOD11_L2 and MYD11_L2 at 1-kilometre resolution obtained at 5-minute satellite collection intervals; other LST products are also available), which has been validated on a number of occasions. It is believed that the MODIS LST archive provides an accurate product with the norm of four images per day (Terra in the morning and evening, Aqua at noon and midnight), thus it was decided to fit SEVIRI LST to MODIS LST in order to improve the SEVIRI LST. Stisen et al. (2004) already fitted SEVIRI to the MODIS LST in order to calibrate their algorithm of retrieving SEVIRI LST in a semi-arid environment in West Africa. They used linear regression with a 5 km MODIS LST product as their reference data. Their regionally based algorithm was tested for three days during the rainy season and the results showed a RMSE of 2.9°C. Before SEVIRI LST is fitted to MODIS LST, the sources of the differences between these products have to be considered. The positive LST difference between SEVIRI and MODIS (Madeira et al., 2005) are most likely to be a consequence of the difference in the zenith angle between the surface and the two sensors. SEVIRI is on board MSG, which is a geostationary satellite situated above the equator, thus it always “sees” Europe from the south. This means, that it sees more of the southern slopes and less shadows behind the northern slopes as it would if it saw the surface from the zenith. On the other hand, MODIS is on board the polar orbiting satellites Terra and Aqua that can “see” the surface at different zenith angles, thus it is able to see more shadows (where the LST is lower than in insulated areas) than SEVIRI does (figure 15). Moreover, this difference can be a consequence of relief shades as well as Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. land cover shades – an insulated part of the land cover is warmer than a non-insulated one, and this could explain the seasonal bias (Noyes et al., 2006). The vegetation has a higher density during the summer, thus SEVIRI cannot see through it and detect the warmer southern parts of the plants during the summer. Finally, the LST generation from MODIS and SEVIRI images is different; the most important difference is using the different emissivity data (LandSAF uses FVC and land cover data in order to compute LST from SEVIRI images, while NASA uses quarterly land cover and daily snow cover to compute LST from MODIS images). Sunlight produces many shadows during the morning overpass of the MODIS, which SEVIRI cannot detect. MODIS on TERRA Equator SEVIRI on MSG Figure 15: SEVIRI always looks upon Europe from the south, thus it sees less shadows (in red) and a greater part of the southern area than MODIS. Slika 15: SEVIRI gleda proti Evropi vedno z juga, zato vidi manj senc (rdeče) in več južnih pobočij kot MODIS. The easiest solution was to fit SEVIRI to MODIS LST by multiple regression, however MODIS LST had to be first of all averaged within a 5 km window (5×5 pixels). Additional explanatory variables can be used beside SEVIRI LST: Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. • solar declination and solar zenith angle, • albedo, • FVC (this is also a LandSAF product but it is still in its demo development status), • elevation, • correction regarding relief morphology and land cover. FVC was used to correct the differences in emissivity and solar declination to correct the seasonal cycle. Furthermore, the difference distribution corresponds to the elevation and albedo distribution. The correction regarding relief and land cover is a consequence of the geometry relationship among SEVIRI, the Sun and the land surface: SEVIRI always sees the same areas that are sometimes more and sometimes less exposed to the Sun. A number of topographic correction techniques are available (Liang, 2004), however they can only be applied to sensors with a higher spatial resolution. Therefore, a new approach is suggested in the study. 3.2.1. Relief morphology correction factor The SEVIRI has a constant azimuth and zenith angle above the surface, thus one can compute its incidence angle, which measures the quantity of the relief the SEVIRI can see, from its zenith angle and azimuth (figure 16). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. MSG ..– .MSG local meridian North Pole standing point MSG .standing point equator equator local parallel local meridian zMSG Figure 16: The position of MSG SEVIRI in the equatorial (left) and meridian (right) plane. Slika 16: Položaj MSG SEVIRI v ekvatorski (levo) in poldnevniški (desno) ravnini. The zenith angle can be computed with the scalar product of the vector that defines the normal to the Earth (represented by the ellipsoid) in the standing point, and the vector that defines the direction of the satellite from the standing point. If both vectors are normalised (their length equals one) and defined within a geocentrical coordinate system, the cosine of the angle between them (zenith angle) equals their scalar product. A similar procedure can be applied to compute the azimuth of the satellite, but both vectors have to be first projected to the equatorial plane. If one assumes that solar radiance exposure has a large influence on LST, the solar incidence angle ., which is the major influence of relief morphology on the solar radiance exposure, has to be considered too. Once more, this can be done by vectors; it is more common to use equations derived from spherical geometry (Allen et al., 2006). .+(t -12) .= ·.+e 12 =sin()sin() +cos . cos() ·cos . cos z .·. ()·. () (3) cos()·sin() . cos h -sin()cos . . ·() .·() cosa = sin() z cos.=cos() ·cos()+sin .·sin()z ·cos a -µ) ß z () ( Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. In both cases, the inclination ß and the aspect µ of the slope (computed from DEM), geographical position (., .) of the standing point, declination ., equation of time e (e.g. Baldocchi, 2006) and daily time t are required as input data. The equation (3) does not consider cast shades (behind an obstacle), thus it is appropriate only for low spatial resolution. At the end one can estimate the spatial distribution of the solar radiance exposure received by the ground seen by SEVIRI by multiplying the solar incidence angle and the SEVIRI incidence angle. 3.2.2. Land cover correction factor A detailed land cover model containing the elevations of the vegetation and anthropogenic objects above the terrain is usually not available, thus it is difficult to estimate the geometry between SEVIRI, the Sun and the land cover surface. The roughness of the land cover surface, which causes micro shadows, has a large influence on LST. LST is too high over rough surfaces, because SEVIRI sees a higher fraction of plants and anthropogenic objects facing south (figure 17), thus the LST of the pine forest has to be corrected while the LST of grassland does not. Figure 17: The zenith angle of the satellite on the surface has a significant influence on LST. If the sensor records an image from the side (left), the image is deformed, thus the sunny areas might be overrepresented compared to the shades, which is not a problem with the images recoded from the top (right). Slika 17: Kot gledanja satelita na površino ima pomemben vpliv na temperaturo površja. Če senzor beleži podobo od strani (levo), je ta popačena, zato so lahko osončena območja močneje zastopana kot sence, do česar ne pride pri navpičnih posnetkih (desno). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The influence of the surface roughness can be estimated by applying an average surface shape to each land cover class. One can assume that land cover surface consists of two parts – a surface facing the satellite sensor and a surface facing the opposite direction. The sensor always sees the surface orientated towards it larger than it actually is. The largest affect on the image occurs when the surface is perpendicular to the line of the satellite view. The parts of the pixel that are covered by the surface looking towards and away from the sensor depend on the inclination of these surfaces and can be written in a mathematical form (see figure 18 to clarify the symbols). ß.MSG ..MSG – .v 1 pixel d 1 – d b a c Figure 18: Influence of the inclined surface upon the satellite image – a sensor detects more of the surface that is orientated towards the sensor (light coloured). Slika 18: Vpliv nagnjene površine na satelitsko podobo – senzor zazna več tistega dela površine, ki je obrnjen proti njemu (svetlo obarvano). v = a · tan.= b · tan ß= c · tan. MSG . .tan ß· tan. 1 .. v = b + c =.cos.MSG ·(tan ß· tan.) cos. MSG . (4) vv tan.. 1 tan ß. d = b + a =+= ·.. +.. tan ß tan. tan ß+ tan. cos. sin. MSG . MSG MSG . Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. If the surface looking towards the Sun and the opposite direction are presented as the edges of triangles (figure 18), the portions (d and 1 – d) can be derived from the height of the triangle v. The solar inclination angle of the surface has to be estimated in order to estimate the geometrical relationship between SEVIRI, the Sun and the land cover surface (figure 19). Two assumptions have been made in order to simplify the calculation: • one part of the surface is orientated towards the SEVIRI (it has the same azimuth as SEVIRI) and the other one in the exact opposite direction a = a +. = a , b c M • the surface is symmetrical – both parts of the surface have the same inclination ß ß =. . ß..S – ..b . c .S Figure 19: Influence of the inclined surface on the solar radiance exposure. Slika 19: Vpliv nagnjene površine na osončenost. The assumption as regards symmetrical objects representing the land cover might not be reasonable at first sight, because nothing in nature is perfectly symmetric. However, an average tree canopy can be represented with a cone or a paraboloid. If one also considers the viewing angle of SEVIRI, it is clear that SEVIRI cannot observe a surface that is parallel to the SEVIRI line of view. Therefore, it makes sense to represent the land cover surface only with objects made of surfaces that are perpendicular to the SEVIRI line of view. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Because of the above assumptions, the incidence angles on both parts of the surfaces can be written as the next equation. cos()=cos() ·cos z +sin()·sin(z)·cos(-aM .d ß() ß aS ) (5) () 1 d cos() ·cos()sin() () · z cos -aM ) cos .-= ß z -ß sin ·(aS Therefore, if LST was influenced merely by solar radiance exposure and the geometry of the land cover, SEVIRI would detect the next signal. L =d ·cos(.)+(1-d )·cos(.d ) (6) Cd 1- Using equation (4), equation (6) can be rewritten as the final expression of the function of the solar zenith angle z and solar azimuth aS; in which A and B are constants that depend merely on the land cover and the elevation angle of the SEVIRI (7). . 1 tan ß. L =. + -1.·sin()ß·sin() ·cos a -a )+cos() ·cos() C.. z (S MSG ß z cos. sin . . MSG MSG . ()·cos(-a )+B ·cos() L =A ·sin za z C S MSG (7) . 1 tan ß. 2d ·sin()-sin()=. + -1.sin ß A = ß ß .·() . .cos.MSG sin .MSG . B cos() =ß A SEVIRI pixel is relatively large, thus it cannot be expected that the surface will be divided merely into two parts within one pixel. Therefore, the surface roughness and the influence of solar radiance exposure on LST have to be first of all estimated in a high spatial resolution (compared to SEVIRI) and then averaged into the SEVIRI spatial resolution. In the end, all correction influences are considered within the multiple regression in order to correct the original SEVIRI LST to fit the MODIS LST, which is believed to be more accurate. Only then is it possible to downscale LST to 1000 m spatial resolution and compare the results with MODIS LST of the same spatial resolution. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 3.3 Downscaling Literature on downscaling is quite rare (except when dealing with panchromatic sharpening – a fusion of high spatial resolution panchromatic images with low spatial resolution colour images), especially in the case of LST. However, the examples of downscaling ambient air temperature (Agnew and Palutikof, 1996; Chung and Yun, 2004; Kostopoulou et al., 2006) showed that numerous methods can be used for downscaling LST: • co-kriging is a geostatistical interpolation technique that uses additional data correlated to the data to be downscaled but it requires a cross-variogram and it is appropriate only for point data, • neural networks provide a powerful tool when sufficient data (sufficient reference and explanatory data) is available, • wavelet transform has already been used for panchromatic sharpening but it can fuse only two similar images, and • contextual analyses (part of the map algebra analyses) can be used to determine the local relationship among the different spatial variables in low spatial resolution and use these relationships on the high spatial resolution data. Various authors (Saravanapavan and Dye, 1995; Prihodko and Goward, 1997; Czajkowski et al., 2004) have documented that LST is linearly correlated to the normalized difference vegetation index (NDVI; figure 20). The relationship is not constant – it varies with time and space. Therefore, the optimal way of establishing the correlation parameters is by applying a linear regression analysis between LST and NDVI within a moving window analysis (a contextual analysis). If the correlation parameters are established from low spatial resolution data, these parameters can be used for LST downscaling by applying these parameters to high spatial resolution NDVI (figure 21). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 20: Influence of NDVI on LST [°C] during the daytime [h] (darker colours present points with a higher NDVI; 29.08.2005) on the chosen stations marked with different symbols. As expected, it can be seen that the areas covered in vegetation do not heat up as much as those with less vegetation. Slika 20: Vpliv NDVI na temperaturo površja [°C] prek dneva [h] (temnejše barve prikazujejo točke z višjim NDVI; 29. 08. 2005) na izbranih postajah, označenih z različnimi simboli. Po pričakovanjih lahko vidimo, da se območja z bogatejšim rastjem ne segrejejo toliko kot manj poraščena območja. The same procedure can also be applied to additional data if the spatial distribution of LST can be explained with them; for example, a visual comparison of the elevation and LST distribution implies that the elevation is also highly correlated to LST. The elevation is a complex parameter due to some special weather situations such as temperature inversion. Therefore, it cannot be said that higher elevation always leads to lower LST, but it can be predicted that the relationship between LST and elevation changes gradually over space, thus a regression made by a contextual analysis should be a suitable tool to evaluate the relationship. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. LST LST = k·NDVI + n 1 regression on low resolution data within the moving window in order to determine the correlation between NDVI & LST NDVI LST 2 insert the high resolution NDVI into the equation NDVI 3 compute the high resolution LST NDVI LST Figure 21: Downscaling with a high spatial resolution NDVI. Slika 21: Izboljšanje ločljivosti z uporabo NDVI visoke prostorske ločljivosti. A visual identification of the correlation with other appropriate data has to be performed too. The correlations between this data and LST are computed for each low resolution pixel within the moving window around it and estimated whether the correlations are not a matter of chance. All significant data is then used within the multiple regression. Each low resolution pixel is assigned regression parameters, which are resampled from the original low spatial resolution data to the resolution of the additional high spatial resolution data. At the end these Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. parameters can be applied to the high spatial resolution data, which yields into the LST spatial resolution improvement. 3.4 Case study: Slovenia 3.4.1. Fit of SEVIRI to MODIST LST This approach was tested in Slovenia. At first fourteen LST MODIS images were chosen to fit the SEVIRI to MODIS LST. The images had to fulfil two conditions: the zenith angle of the MODIS had to be significantly smaller than the zenith angle of the SEVIRI (it makes no sense to make a topographical correction and use images with topographical distortions as a reference) and at least one third of the area had to be cloudless. Not many appropriate images meeting these conditions were available for the 2005–2006 period (the LandSAF data is available since March 2006). All the appropriate images from Terra and Aqua done by day and night were used in order to cover different seasons and different daily times (table 6). Table 6: : Images used to fit the SEVIRI LST to MODIS LST; there are five images with a zenith angle greater than 25°, which is approximately half of the SEVIRI zenith angle. Preglednica 6: Podobe uporabljene za prilagoditev SEVIRIjeve na MODIsovo temperature površja; pet podob ima zenitni kot večji od 25°, kar je približno polovica zenintega kota SEVIRIja. Day of the year (date) Day/night image (GMT) Satellite Zenith angle [°] Shift in x (+ due E) Shift in y (+ due S) 258 (15.09.2005) night (22:15) Terra 33 3 2 302 (29.10.2005) day (12:15) Aqua 25 2 3 313 (09.11.2005) day (12:00) Aqua 9 5 -1 338 (04.12.2005) day (10:00) Terra 14 0 -5 044 (13.02.2006) night (00:45) Aqua 30 0 -3 097 (07.04.2006) day (12:15) Aqua 25 2 4 131 (11.05.2006) day (12:00) Aqua 8 3 3 137 (17.05.2006) day (09:45) Terra 31 -2 -5 167 (16.06.2006) day (10:00) Terra 8 -1 -3 200 (19.07.2006) day (12:15) Aqua 34 3 2 236 (24.08.2006) night (01:00) Aqua 28 -1 -1 245 (02.09.2006) day (12:00) Aqua 12 4 1 291 (18.10.2006) night (01:00) Aqua 21 -2 -2 323 (19.11.2006) night (01:00) Aqua 15 0 0 SEVIRI and MODIS images were also additionally geo-referenced before the fit began, because of the systematic shift of SEVIRI images (they all seemed to be moved one pixel too far to the North) and the unsystematic shift of the MODIS images (shift parameters in pixel Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. size of one kilometre are listed in table 6), which were observed by comparing LST with the higher spatial resolution data. 3.4.1.1 Relief morphology correction factor In order to consider the effect of relief morphology a constant zenith angle and azimuth of SEVIRI were computed for entire Slovenia (the variable zenith angle and azimuth do not significantly change the results, since the area is rather small and the relief is rough). It was computed that the zenith angle of the MSG equals 49° and the azimuth 200° – these values were used to compute the SEVIRI incidence angle in Slovenia (figure 22). Figure 22: Incidence angle of SEVIRI in Slovenia averaged to 1 km pixel (elevation data source: DMR 100, November 2005, © Surveying and Mapping Authority of the Republic of Slovenia). The cosine of the angle is a measurement of the surface exposure to SEVIRI. SEVIRI sees more of those areas that are more exposed to it (the incidence angle equals between 40° to 50° for most of the area). Slika 22: Vpadni kot SEVIRI-ja v Sloveniji povprečen na slikovni element velikosti 1 km (vir višinskih podatkov: DMR 100, november 2005, © Geodetska uprava Republike Slovenije). Ta kot je mera za izpostavljenost površja SEVIRI-ju. SEVIRI vidi več tistih območij, ki so mu bolj izpostavljena (vpadni kot znaša 40° do 50° na večini območja). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The source of the relief elevation was the digital elevation model DMR 100 produced by the weighted sum from various sources (Podobnikar and Mlinar, 2006). The incidence angle of SEVIRI, which is the stationary part of the relief morphology correction, was averaged to 1000 m spatial resolution by a median within a 9 by 9 moving window. The solar incidence angle computed by equation (3) was taken into consideration by multiplying it with the incidence angle of SEVIRI in order to complete the relief morphology correction. The correction was preformed at 1000 m spatial resolution and then averaged to 5000 m spatial resolution by a median within a 5 by 5 moving window. 3.4.1.2 Land cover correction factor The Corine land cover 2000 data (EEA, 2005; it has poorer spatial resolution than the data used in chapter 2, which does not cover the entire study area as the Corine data does) was reclassified according to the roughness of the land cover surface (described in chapter 3.2.2; figure 23). The values of the average inclination angle applied to each land cover class were estimated intuitively: the coniferous forest and bare rocks have the largest average inclination (50°), natural grasslands, water bodies, etc. are flat (0°) and the other land cover types have the average inclination somewhere in between these values. Two datasets were derived, indicated as A and B because they correspond to parameters A and B from the equation (7). Both layers were averaged to 1000 m spatial resolution by a median within a 9 by 9 moving window. The correction was completed using the solar zenith angle and azimuth calculated with equation (7) at 1000 m spatial resolution and then averaged to 5000 m spatial resolution with a median within a 5 by 5 window. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 23: Layers A (above) and B (below) used for land cover correction (data source © EEA, Copenhagen, 2005); look at equation (7). Slika 23: Sloja A (zgoraj) in B (spodaj) za popravek pokrovnosti (vir © EEA, Copenhagen, 2005); glej (7). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 3.4.1.3 Calibrating the SEVIRI LST correction It was observed that the average SEVIRI LST ( LSTSEVIRI ), the cosine of the solar zenith angle zs, the solar declination ., the relief morphology correction RC, the land cover correction LC, elevation EL, albedo AL and the fraction of vegetation cover FVC are correlated to the LST differences between SEVIRI and MODIS. The last three have correlation indexes (with MODIS LST) smaller than 0.2 (all others larger than 0.7), which seems statistically insignificant, but combining their correlations with other parameters improves the results (the standard deviation of the fit improves by 15%). Therefore, the SEVIRI LST ( LSTSEVIRI ) was fit to the MODIS LST through the multiple regression from the listed parameters. The images from table 6 were used to calibrate the general correction for the SEVIRI LST. LST =-1.47 + 0.688· LST + 0.285· LST - 6.48· cos(z )- SEVIRI SEVIRI S - 5.75· . - 2.20 ·10-3 · h + 8.85· Rc - 0.82 · Lc -1.43· AL + (8) -4 -72 + 8.89 ·10 · FVC -1.12 ·10 · FVC Figure 24 shows the scatter of SEVIRI LST regarding MODIS LST averaged to 5000 m spatial resolution; the correlation is relatively high (r = 0.98) but the extreme differences can be up to 15 °C. Figure 25 shows the same scatterplot after the correction of SEVIRI images. The comparison between figure 24 and figure 25 shows that the proposed correction improved the LST model – the correlation is higher (r = 0.99), the standard deviation from the model equals 1.6 °C and the extreme differences is significantly smaller (min = -9.7 °C; max = 8.9 °C), although they are still relatively large. RMSE varies for each day but in general it is around 1.5 °C (between 0.9 °C the 29.10.2005 an 2.8 °C for the 17.05.2005). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 24: Scatterplot MODIS-SEVIRI LST for the days listed in the table 6 (data for each day is in a different colour). Slika 24: Razsevni grafikon med MODIS-ovo in SEVIRI-jevo temperaturo površja za dni iz preglednice 6 (vsak dan je na grafikonu drugače obarvan). Figure 25: Scatterplot MODIS-corrected SEVIRI for the days listed in the table 6 (data for each day is in a different colour). Slika 25: Razsevni grafikon med MODIS-ovo in popravljeno SEVIRI-jevo temperaturo površja za dni iz preglednice 6 (vsak dan je na grafikonu drugače obarvan) v letih 2005 in 2006. Figure 26 presents an example of the differences between the corrected SEVIRI and MODIS LST. 60% of all the differences are smaller than 1 °C in the absolute sense. The images for the other days do not show a systematic spatial distribution of the differences, thus the rest of the variation can be a consequence of random influences, such as different times of recording Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. the LST (the SEVIRI makes an image every 15 min, thus the time of recording can differentiate by a maximum of 7.5 min). The only systematic difference that was not removed with the shown procedure is the different “contrast” within MODIS and SEVIRI LST. It seems that MODIS has a greater range within its LST, thus the inclination of the fitted line through both data is different – the line fitted through SEVIRI data is never as steep and this drawback is not solved with the proposed procedure – the inclination differences of the fitted lines remain after the correction. Figure 26: Difference between corrected SEVIRI and MODIS LST; an example for 07.04.2006 at 12:15 GMT. Slika 26: Razlika med popravljenimi SEVIRI-jevimi in MODI-sovimi temperaturami površja; primer za 07. 04. 2006 ob 12:15 GMT. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 27: A scatterplot MODIS-SEVIRI LST before the fit to MODIS LST; an example for 29.10.2005 at 12:15 GMT. Slika 27: Razsevni grafikon med MODIS-ovo in SEVIRI-jevo temperaturo površja pred prilagoditvijo na MODIS-ove temperature; primer za 29. 10. 2005 ob 12:15 GMT. Figure 28: A scatterplot MODIS-corrected SEVIRI after the fit to MODIS LST; an example for 29.10.2005 at 12:15 GMT. Slika 28: Razsevni grafikon med MODIS-ovo in SEVIRI-jevo temperaturo površja po prilagoditvi na MODIS-ove temperature; primer za 29. 10. 2005 ob 12:15 GMT. 3.4.2. Downscaling from 5000 to 1000 m spatial resolution The visual comparison of MODIS LST and some of the available data suggested that there are many explanatory variables (of 1000 m or even higher spatial resolution) available for Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. downscaling. The MODIS vegetation indices product (MOD13A2) does not contain merely NDVI but also the enhanced vegetation index (EVI), which has a higher contrast in the highly vegetated areas and can be used also for downscaling. The artificial vegetation index was generated from the reclassified Corine 2000 land cover data (EEA, 2005; 0 for artificial surfaces and 1 for coniferous forest, other classes in between) in order to provide a higher contrast during certain periods when NDVI and EVI can be similar all across the area. Moreover, the LST spatial distribution corresponds well to the albedo, which is another MODIS product (16-day composite; MOD43B). As expected, LST also corresponds to the elevation and solar incidence angle; both values were derived from the digital elevation model DMR 100 (Podobnikar and Mlinar, 2006). The correlation between LST and NDVI was proven already in the past (e.g. Prihodko and Goward, 1997). Therefore, it was decided to always consider NDVI in the downscaling procedure. All other explanatory variables are considered only if their correlation is statistically significant (not a consequence of chance). The correlation between LST and the explanatory variables was always established within a 7 by 7 moving window (the correlation was also tested in a 5 by 5 window, but the results were not as good) and in the event that at least 16 pixels (one third of the moving window) with retrieved LST were available, the correlation index (between LST and each explanatory variable) was computed and tested according to equation (1). The significance interval of 0.9 was used, which means that the correlation index at 14 degrees of freedom should always be larger than 0.34 in order to consider the explanatory variable within the further analysis. Figure 29 presents an example of corrected SEVIRI LST. The image was used as input data in the downscaling procedure (an example of a validation and the resulting LST are in the figure 30). Note that the no-data areas enlarge during the downscaling process because of the characteristics of the moving window analysis. The correlation index of the results (compared to MODIS LST) is in general greater than 0.8, RMSE varies for each day but it is usually around 2.0 °C (between 1.4 °C the 29.10.2005 an 3.0 °C for the 17.05.2005). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 29: Corrected SEVIRI LST (5000 m spatial resolution); an example for 19.07.2006 at 12:15 GMT. Slika 29: Popravljene temperature površja zaznane s SEVIRI-jem; primer za 19. 07. 2006 ob 12:15 GMT. A remark as regards the validation procedure has to be made because it was not preformed with completely independent data. Only 14 LST datasets were used (table 6) in order to fit SEVIRI to MODIS LST, thus no data was used merely for the validation of the SEVIRI correction. Furthermore, the same 14 MODIS LST datasets were used for validating the downscaled results, because (as already stated) no other data was available. The downscaled LST is partially connected to the data used for the validation since the input for the downscaling is a corrected SEVIRI LST that contains some information also from the 14 MODIS LST datasets averaged to 5000 m spatial resolution. This means that the validation is not performed on totally independent data, but the effect of the explanatory variables has a significant impact on the LST during the downscaling procedure, thus the effect of the MODIS LST, which slightly changes the input low resolution LST, is not to be considered significant. Therefore, the average RMSE of 2.0 °C is an appropriate measurement of the result quality until further data is made available. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 30: MODIS (at 12:20 GMT) and downscaled SEVIRI LST (at 12:15 GMT) for 19.07.2006. Slika 30: Z MODIS-om zaznana (ob 12:20) in modelirana (ob 12:15) temperatura površja za 19. 07. 2006. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The problem that was unsolved by the SEVIRI LST correction (different inclination of the fitted lines) remains unsolved also after the downscaling procedure. In fact, the inclination difference got even higher (because of the low contrast within the original data), thus additional correction was applied to the downscaled LST – it was estimated that the average fitted line (from the downscaled data) has the inclination coefficient of 0.6 instead of 1, thus the downscaled data was additionally rotated around the mean data for each day. Figure 31 shows an example of scattering between MODIS and the downscaled SEVIRI LST; the inclination difference still exists, however, it is no longer systematically smaller than the fitted line inclination of the MODIS LST. Figure 31: A scatterplot MODIS-downscaled SEVIRI; an example for 07.04.2006. Slika 31: Razsevni grafikon med MODIS-ovo in izboljšano SEVIRI-jevo temperaturo površja; primer za 07. 04. 2006. 3.5 Discussion The described downscaling procedure is a simple yet efficient method. It allows one to use numerous different datasets that influence the LST spatial distribution. It is also rather simple to use additional data– it is more problematic to find suitable data that might explain the LST spatial distribution. The method works as expected well when the SEVIRI records LST for areas that are not similar – if the areas are too similar, no contrast within the input data is Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. seen, and this yields into an unreliable regression analysis. Therefore, one cannot easily evaluate the results, because LST is not equally explained with the explanatory variables all across the working area. It can be seen that the areas, in which LST is more generalised – less detail compared to MODIS – correspond to the areas with a low multiple regression correlation index. Therefore, some other explanatory variables should be integrated into the procedure, for example geological data, but unfortunately they were not available for the study. Another disadvantage of the method is its relatively low speed. Moreover, the reference data (MODIS) makes the evaluation of the results even more uncertain. One would expect that the data is accurately positioned (at least data retrieved by MODIS on Terra, which has an advanced navigation system) but they were not. The MODIS data also depend on the viewing angle of the sensor, thus they can not be considered as an error free reference data. When one considers that time of retrieval can differentiate for a couple of minutes and that sensors have different characteristics, then it is a difficult task just to compare the MODIS and SEVIRI data already at 5000 m resolution. The differences at 5000 m spatial resolution are also the consequence of different emissivity data; the emissivity differences were not available for the study, thus FVC was used as a their proxy. The described approach of fitting SEVIRI to MODIS LST is promising although it has to be tested in another case study area because the used parameters are most likely to be representative only for the presented case study area. Furthermore, the MODIS LST datasets used in the study is not a completely independent validation data because the same datasets were used already in the fit of SEVIRI to MODIS LST. An alternative validation option would be point measurements of LST with spectrometers or airborne spectrometer measurements. However, this data was not available in the study. On the other hand, measurement of the temperature under the surface (depth of 2, 5, 10, 20 and 30 cm) were available on some meteorological stations, but the differences between LST and ground measurements can be significant because of a high temperature gradient just under the surface. Moreover, satellite usually does not observe bare ground but ground covered by vegetation and anthropogenic objects. The research in the LST downscaling field has been rare so far, thus the results cannot be truly compared to similar studies. Nevertheless, the described downscaling procedure Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. provides LST in 1000 m spatial resolution every 15 min with an average RMSE of 2 °C. This is a great improvement considering either the poor spatial resolution of SEVIRI (5000 m with the same accuracy) or the poor temporal resolution of MODIS (combining Terra and Aqua day and night images 6 h in average; but with an estimated accuracy of 1 °C). In the near future, the method will also be tested with independent data (when appropriate data will be available). If the RMSE of 2 °C according to MODIS LST remains, than one can expect an absolute accuracy (according to the error propagation; Heuvelink, 1998) of 2.2 °C. One can consider the resulting LST as a good funding to parameterize ambient air temperature from it (that is the main purpose of the thesis), but its spatial resolution might still be too coarse for high mountains. Therefore, additional downscaling procedure is described in the next chapter. Many different influences on LST were considered within this chapter, but in the next one it is assumed that solar radiance exposure is the driving factor for the differences in LST. LST spatial resolution is then improved by fusion of LST and solar radiance exposure of eight times better spatial resolution by wavelet transform. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. MODELLING LAND SURFACE TEMPERATURES WITH WAVELET TRANSFORM One of the downscaling possibilities mentioned in the previous chapter is the wavelet transform, which transforms the original data into a set of coefficients that are often easier to process than the original data. Some examples of image fusion based on wavelet transform are available in the remote sensing field. With the integration of spatial and spectral information from different sensors, multi-sensor image fusion provides a composite image with a more complete and accurate description of the scene than any of the original source images (Ranchin et al., 2003; Švab and Oštir, 2006). Wavelet transform has proven to be an appropriate approach also for downscaling since it can improve the spatial resolution of low resolution data with a minimal distortion of the original image’s spectral content (Acerbi- Junior et al., 2006). The wavelet transform has been mainly applied to improving the spatial resolution of a multi-spectral image by merging it with a panchromatic image. Some examples include the fusion of images acquired with different sensors (Landsat TM 5 and SPOT); some include the fusion of images acquired through different channels of the same sensor (Ikonos; Ranchin et al., 2003). The ratio of the spatial resolution between both source images is usually 2 to 4, as can be seen in the case of the Ikonos multispectral image (4 m) and panchromatic image (1 m). However, the fusion of source images with a larger ratio of spatial resolution is also possible – MODIS (250 m and 500 m spatial resolution in the Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. channels 1–7) and Landsat TM 5 have been successfully fused by wavelets transform, which means that even the ratio of 16 is possible (Acerbi-Junior et al., 2006). All of these examples are based on fusing the “related” images – images with similar spectral information. However, the purpose of this study is to downscale the land surface temperature (LST) with the data that influences it and is available at a higher spatial resolution. This means that the images with completely different information have to be fused, thus it is important to understand that LST varies with land and relief type – it seems that the solar radiance exposure represents the greatest influence in the case of a rough relief while soil humidity represents the greatest influence in the case of a mild relief. Therefore, the optimal option is to model LST in high spatial resolution using data on specific thermal capacity of the material, wind speed, solar insulation, etc. However, the correlation between these parameters has not been investigated in details, thus it is difficult to estimate the LST. As Slovenia is a rough country, it was decided to investigate the influence of the solar radiance exposure on LST; in general the south facing slopes receive up to 10% more and the north facing slopes up to 80% less solar irradiance as a horizontal surface on a yearly level (Kastelec et al., 2007). It can be seen (especially in spring or autumn mornings) that incoming solar radiance exposure has a significant influence on LST – areas in the relief or vegetation shadow often tend to be covered by frost, while areas exposed to the Sun tend to be warmer, thus they are not covered by frost (figure 32). If the weather is cloudy (the sky is covered by clouds during approximately 50% of daytime in Slovenia), other parameters are more important. Furthermore, the correlation between solar radiance exposure and ambient air temperature was examined in a previous study (Chung and Yun, 2004); since air is mainly heated indirectly by the surface, it can be assumed that the incoming solar radiance exposure is the parameter that has the largest influence on LST in mountainous areas with homogeneous land cover. Therefore, simulated solar radiance exposure data was fused with the LST obtained by downscaling SEVIRI LST to 1000 m spatial resolution in order to downscale LST to an even higher spatial resolution – 125 m. The downscaling was preformed in two steps for validation purposes, which can be controlled also by MODIS LST in 1000 m spatial resolution. Moreover, the first downscaling step makes it possible to consider a greater Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. number of influences, while the second step takes in consideration only solar radiance exposure. Figure 32: Direct solar radiance exposure is much stronger than the diffuse radiation, thus the areas in the shade have a significantly lower LST (frost is only in the shade). Slika 32: Direktna osončenost je mnogo močnejša od difuzne, zato imajo območja v senci nižjo temperaturo površja (slana je le v senci). 4.1 Discrete wavelet transform In order to build representations of signals that are able to adapt to the nature of the signal, scientists over the last century tried to overcome the limitations of the Fourier transform, which is probably the most popular transform for signal processing. Further information on Fourier transform and continuous wavelet transform can be found in Appendix B, thus merely a brief outline of the discrete wavelet transform is found below. Wavelet transform is a mathematical transform that is applied to signals in order to obtain further information from the signal, information that might not be immediately visible from the raw signal itself. It decomposes a signal through a family of wavelet functions obtained from the translation and dilatation of the wavelet function (Pajares and Manuel de la Cruz, 2004). The continuous wavelet transform requires an abundance of computing resources, thus the discrete wavelet transform was developed, in which the wavelets are discretely sampled. The discrete wavelet transform also makes it possible to combine certain coefficients deriving from various images (Pajares, Manuel de la Cruz, 2004). The discrete wavelet transform uses filters with different cut-off frequencies to analyze signal x at different scales: in order to Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. analyze the high and low frequencies the signal is passed through a series of high and low- pass filters. First, signal x is passed through a low-pass filter with impulse response g. Filtering the signal corresponds to the mathematical operation of convoluting the signal with the filter impulse response. The convolution operation y in discrete time is defined with the following equation. . y[]n=.x[] k · g[n - k] (9) k =-. Signal x is also decomposed using a high-pass filter h. The output provides the detail and approximation coefficients (from the high and the low-pass filter). It is important that the two filters are related to each other. A half band low-pass filter removes all frequencies that are above half of the highest frequency within the signal. Since half of the signal frequencies are removed, half of the samples can be discarded. The filtered outputs are then downsampled by 50%, simply by discarding every other sample. This constitutes one level of decomposition and can be expressed mathematically with the following equation (10). low []=.. [ ] · g 2· n - k] yn xk [ k =-. yhigh []n =.[] [ · n - k . xk · h 2 ] (10) k =-. This decomposition halves the time resolution because only half of the samples describe the entire signal (Pajares and Manuel de la Cruz, 2004). On the other hand, it doubles the frequency resolution, since the signal’s frequency band spans only across half of its previous band. As seen below on the decomposition-reconstruction diagram the procedure can be repeated for further decomposition until level 3. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. x[n] g[n] h[n] ˇ2 ˇ2 Level 1 coefficients g[n] h[n] ˇ2 ˇ2 Level 2 coefficients g[n] h[n] ˇ2 ˇ2 Level 3 coefficients x[n] g -[n] h-[n] ^2 ^2 Level 3 coefficients g -[n] h-[n] ^2 ^2Level 2 coefficients g -[n] h-[n] ^2 ^2Level 1 coefficient Filter diagram of decomposition and reconstruction until level 3: ˇ2 (^2) stands for the downsampling (upsampling) operator by the factor of 2 (Chen and Zhu, 2007). Diagram filtra dekompozicije in rekonstrukcije do 3. stopnje: ˇ2 (^2) predstavlja operator za zmanjšanje (večanje) podrobnosti vzorca za faktor 2 (Chen in Zhu, 2007). The reconstruction follows the decomposition procedure, just the other way around. The signals are upsampled by two on every level, passed through the low-pass synthesis and high- pass filters and then added. The analysis and synthesis filters are identical, except for the time reversal. [] . y [] [ ·k -n]+y [] [ k ·h 2 ·k - xn =. k ·h 2n] (11) low high k =-. The described procedure is valid for a one dimensional signal. A two dimensional image can be treated as a group of one dimensional signals – in this case the procedure is applied to all rows, all columns and all diagonals. In order to fuse the two images, another consideration has to be taken into account after both images are transformed – i.e. one has to know how to merge wavelet coefficients. A weighted averaged sum is always a good option, but when images with different input resolutions are used, one can also merely exchange some of the coefficients from the high resolution image with the corresponding coefficients from the low resolution image (Pajares and Manuel de la Cruz, 2004). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 4.2 Solar radiance exposure simulation Solar energy drives the Earth’s climate and even slight variations in solar radiance could offset or increase global warming. The 23.5° angle between the Earth's spin axis and its orbit around the Sun creates the seasonal cycle, causing the day length and the sunlight angle to vary during the year. As a result, summer is much warmer than winter, and the polar regions are significantly colder than the tropics. The Earth orbits around the Sun on an elliptical path, at which the Earth is slightly closer in early January than it is in July. A consequence of this is that the Earth receives slightly more solar irradiation in January compared to July. Both seasonal effects, the tilt of the Earth’s axis, and the orbital distance to the Sun, are stable and predictable. The solar irradiation received by the surface is influenced by (e.g. Zakšek et al., 2005): • Sun position regarding the Earth, • morphology of the Earth surface, and • climate and transmisivity of the atmosphere. The geometry between the Sun and the Earth is predictable, thus it is easy to model its influences if the exact time is known (the hour of the day can be used to compute the solar azimuth and zenith angle when the day of the year is known – this makes it possible to compute the solar declination). The relief morphology influences the incoming solar radiance exposure (figure 33) predominantly by the incidence angle of the Sun, which is determined by the normal to the surface tangent plane. This depends on the relative positions of the Earth and Sun, the geographical position of the point on the surface, slope and exposition in the chosen point. The incidence angle is easy to estimate if the normal vector to the surface and the vector of the sunrays have been previously determined. The relief morphology influences the effective duration of direct solar radiance exposure, which is the period of sunshine in clear weather, subtracted by the period of shade due to relief obstacles. The surface exposed to the Sun does not only receive diffuse radiation, like the surface in the shadow, but also the energy of direct radiation. It is therefore important to know which part of the surface is in the shadow (Zakšek Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. et al., 2005). Furthermore, the relief morphology is also important for the diffuse part of the radiation – the diffuse radiation is reduced by the sky-view factor, which is the visible proportion of the bright sky that contributes an important part to the solar energy received by the surface (Zakšek, 2006). n s tangent to relief ellipsoid z . horizon line restricted with relief . < 1 . = 1 ideal horizon line Figure 33: The influence of relief morphology on direct (above) and diffuse (below) solar radiance exposure. Slika 33: Vpliv oblike reliefa na direktno (zgoraj) in difuzno (spodaj) osončenost. Clouds and fog can be understood as a filter that increases the diffusion of the incoming radiation throughout the atmosphere. Therefore, the number of cloudy or foggy days should be considered as one of the most significant climate information. The climate is studied by parameters obtained through measurements at meteorological stations. The duration of solar radiance exposure – which is the most important data used for the meteorological model – is commonly recorded with the Campbell-Stokes heliograph. A glass globe focuses the radiation beam onto a special recording paper upon which the trace is burned as the Sun moves across the sky. When the sunrays do not make it through the atmosphere nothing is recorded on the Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. paper. Climate data is important for general prediction (for a longer period). On the other hand, weather data is important in the analysis of a specific period. According to the hypothesis of this study, LST is influenced by the effect of the relief morphology on the solar radiance exposure in mountainous areas. Therefore, only the influence of the relief morphology on solar radiance exposure is modelled and then the influence is considered as a proxy for LST. 4.2.1. Influence of relief morphology on direct solar radiance exposure The Sun’s incidence angle can be calculated from the normal vector to the surface and the vector in the direction of the Sun (Zakšek et al., 2005). Both vectors are defined in the global geocentric Cartesian coordinate system. The input data enables the calculation of the geographical coordinates, ellipsoidal heights, slopes and expositions. The normalized normal vector to the surface is given by the following equation. .- cos.- sin . 0..sin. 0 - cos.. .sin.·cos µ. vn . .. ... =- sin . cos. 0 · 01 0 · sin.·sin µ (12) . .. ... . 0 01..cos. 0 sin... cos.. . .. ... v n is the normal vector to the surface, . is the geographical (ellipsoidal) latitude, . is the geographical (ellipsoidal) longitude, . is the relief slope, and µ the relief aspect. The sunray vectors can be easily defined if the Sun motion is simulated. An observer on the Earth can notice that the Sun's position changes all the time. Therefore, the sunray vector is a function of time and the Sun’s declination. .cos.·sin.. vs .. = cos.· cos. (13) .. . sin.. .. vs is the sunray vector, . is the Sun declination and . is the time angle defined by the following equation (t is time; GMT in hours). t ( = 12 h ) .· - 6 . (14) h Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. One can compute the angle between the normal vector to the surface and the sunray vector by using the scalar product. If both vectors are normalized, then the value of the scalar product equals the cosine of the angle between the selected vectors. If the surface is not exposed to the Sun, it is in the shadow. There are two types of shadow – hill shade and cast shade. The hill shade occurs when the surface is oriented away from the light source. The hill shade can be easily determined if the incidence angle is previously calculated. In contrast, the determination of the cast shade is a more complex procedure. The cast shade occurs only on surfaces that are orientated towards the Sun when there is a topographic obstacle between the Sun and the surface (figure 34). exposed to the Sun partially insolated hill shade cast shade n1 n3 n4 n2 vS vS vS .1 .2 .3 .4vS Figure 34: Solar incidence angle determines whether the surface is in the hill shade. Slika 34: Vpadni kot sonca odloča, ali je površje v lastni senci. The surface could be exposed to the Sun when the incidence angle is smaller than 90° (.1). If the sunrays are almost parallel to the surface (the incidence angle is close to 90°), then the surface is only partially exposed to the Sun due to vegetation, relief micro-morphology and man-made obstacles (.2). If the incidence angle is greater than 90°, the surface is in the shadow – hill shade (.3). Due to the topographic obstacles, some parts of the surface could be in the shadow even if the incidence angle is smaller than 90° – cast shade (.4). Solar radiance exposure can be divided into direct and diffused. The surfaces in the shadows receive only the energy from the diffused radiation, which is often at least twice less than the energy obtained from the direct solar radiance exposure. The first step in determining the cast shade is to find the top of the obstacle. The obstacle's highest point is the point at which the hill shade appears first. This is where the incidence angle becomes smaller than 90°. After finding a potential obstacle (figure 35) the algorithm moves through the digital elevation Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. model (DEM) in the direction closest to the azimuth of the sunrays. The azimuths to the adjacent cells are calculated and compared to the azimuth of the sunrays. –vS vS n0 z0 vS n0 z0 vS n0 z0 vS n0 z0 z2z1 z3 cast shade border exposed to the Sun Figure 35: Cast shade determination. Slika 35: Določanje vržene sence. If the zenith distance from a standing point to the top of the obstacle (z1) is smaller than the zenith distance to the Sun, then the surface is in the cast shade. If the zenith distance to the top of the obstacle (z3) is greater than the zenith distance to the Sun, then the surface is exposed to the Sun. If the standing point is on the border between the illuminated area and the shadow, then the zenith distance to the top of the obstacle (z2) equals the zenith distance to the Sun. Thus, the key to determining the cast shade is in finding the border point. .. h0 - h . vv - arctan. .= arccos(s ·nE ) (15) 2 . d . h0 the obstacle height, h is the standing point height, d is the horizontal distance between the vv top and the standing point, s is the sunray vector and nE is the normal vector to the ellipsoid. r The effect of the relief morphology Rdir on the direct solar radiance exposure Rdir can be computed once all cast shadows are determined and all affected pixels are excluded. Rr = R · cos. (16) dir dir Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 4.2.2. Influence of the relief morphology on the diffuse solar radiance exposure The size of the visible sky is important for diffuse radiation – a plane or a peak receives much more diffuse radiation because it comes from the sky and there is much more visible sky as is in the valley (Zakšek, 2006). The most convenient measure for expressing the size of the visible sky is a solid angle. The solid angle . is proportional to the surface area (S) of the projection of the object onto the sphere centred at that point, divided by the square of the sphere's radius, R. S .=k · (17) R2 The solid angle is also called sky-view factor when the proportionality constant k value equals 1 / 2.. The sky-view factor, the value of which is always between 0 and 1, is easily computed when the visible sky is limited with a cone surface (A – angle between the cone and the horizontal plane). . 2. 2 .= d.·d.=1-sin()A (18) .. 0 A However, perfect geometrical bodies are rare in nature, thus the actual surface is approximated with a cone – the inclination of the cone is approximated with the average elevation angle of the horizon (figure 36). Therefore, the horizon elevation angle has to be determined in multiple directions and its mean value is then used for sky-view factor computation. Not all the diffuse radiation is coming from the sky in the rough relief; some is also reflected from the ground of the opposing slopes. Therefore, the ground albedo has to be considered in that part of the hemisphere, where the sky is not visible – ground-view factor (1 – .). In order r to obtain the effect of relief morphology Rdif on the diffused radiation Rdif the sky-view factor is multiplied (assuming that the diffuse radiation is isotropic) with the diffuse solar radiance exposure and added to the incoming solar radiance exposure (direct and diffuse) multiplied with albedo AL and ground-view factor. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. r R = R ·. +(R + R )· AL ·(1-.) (19) difdif dir dif .A Figure 36: The sky-view factor of the cone depends on the angle between the cone and the horizontal plane. Slika 36: Delež vidnega neba znotraj stožca je odvisen od kota med vodoravno ploskvijo in plaščem stožca. 4.3 Case study: mountainous area on the border between Slovenia and Austria 4.3.1. Reference LST data The case study area is situated in the Alps on the border between Slovenia and Austria. It covers an area of 32 by 32 km. This area was chosen because of the reference data available for validation – a Landsat TM 5 image was available for the central part of Slovenia for 30.07.2005 at 9.30 GMT. The chosen study area was not covered by clouds at the time the image was made. Landsat TM 5 has only one thermal infrared channel with a 120 m spatial resolution. Due to the multiple channels in the thermal infrared spectrum that allows the use of other LST determination methods, and the better spatial resolution (90 m), it would be better to use the ASTER data for reference, but unfortunately no appropriate data was available for Slovenia (covered by clouds, only partially in the area, etc.). The radiance of the sixth TM channel was used to compute LST with the following equation. K2 LST= (20) . K1 . ln.+1. . W . Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. W is the spectral radiance, K1 is the first calibration constant (607.76 Wm-2SI-1µm-1) and K2 is the second calibration constant (1260.56 K; Oštir, 2006). The low resolution LST data was computed as described in the previous chapter for entire Slovenia. The data was re-sampled to 125 m spatial resolution, which corresponds to the fourth level of LST wavelets decomposition in 1000 m spatial resolution (discrete wavelet transform halves the resolution – three times in the present case: 1000 m .. 500 m .. 250 m .. 125 m). Therefore, the Landsat LST was also re-sampled to the resolution of 125 m. 4.3.2. Solar radiance exposure proxy The relief data, which influences the solar radiance exposure, were also resampled to the same resolution as the Landsat LST (figure 37). The source of the relief parameters was DMR 100 (Podobnikar and Mlinar, 2006). Elevation, relief slope and aspect were important for the direct radiance exposure, while the sky view factor, determined in the past study (Zakšek, 2006), was important to simulate diffuse radiance exposure. The diffused radiance exposure also depends on the surface albedo, thus the land cover data (Kokalj and Oštir, 2006) was reclassified according to the MODIS albedo product and then generalized to 125 m spatial resolution. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 37: Relief parameters that influence solar radiance exposure (Alpine area of 60 by 60 km). Slika 37: Parametri reliefa, ki vplivajo na osončenost (območje Alp velikosti 60 krat 60 km). A proxy for the entire solar radiance exposure received by the surface can be determined after the effect of the relief morphology on the direct and diffuse radiance exposure is estimated. However, this can only be performed if the ratio between both parts of the solar radiance exposure is known. Numerous methods can be used in order to determine the approximate ratio . between the diffuse and direct solar radiance exposure (Allen et al., 2006), but within the climatic study data that we used the diffused and direct solar radiance exposure energy measurements on synoptic stations were interpolated into a regular grid with a spatial resolution of 1000 m (Kastelec et al., 2007). The climatic data has an approximately 10 day temporal resolution, thus it was considered appropriate for the estimation of the ratio between the diffused and the direct solar radiance exposure without considering the relief effects (figure 38). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 38: The ratio between the diffuse and the direct solar radiance exposure for the 21st decade (a ten day time period) in the year; it can be seen that the ratio is greater within the high mountains due to convective cloudiness typical for clear summer afternoons – the ratio of the lowland (60%) was therefore used because the sky was cloudless in the study area (Alpine area of 120 by 120 km). Slika 38: Razmerje med difuzno in direktno osončenostjo za 21. dekado v letu; razmerje je večje v visokogorju zaradi konvektivne oblačnosti, ki se pojavi na jasne poletne popoldneve – zato je bilo uporabljeno razmerje za nižine (60 %), kajti nebo na testnem območju ni bilo prekrito z oblaki (območje Alp velikosti 120 krat 120 km). A solar radiance exposure LST proxy was determined by considering the proportional parts of direct and downward diffuse solar radiance exposure from the sky. If . (the ratio between the diffuse and the direct radiation) is known, then the solar radiance exposure LST proxy Rr can be determined from (16) and (19). R .= dif Rdir rr r R =Rdir +Rdif (21) r (.,., ., AL, R )=( .+.·.+(+.·AL ·(-. R =f cos() 1 ) 1 ))·R dir dir r R =(cos().+.·.+(1+.)·AL ·(1-.)) proxy The solar radiance exposure LST proxy Rr is generally a value between zero and two that proxy is linearly correlated to the solar radiance exposure Rr that is influenced by the relief. In order to fuse the solar radiance exposure LST proxy with LST, the radiation proxy histogram (figure 39) had to be matched to the 1 km LST histogram. This was performed through linear scaling (Acerbi-Junior et al., 2006), which yields both data into a similar range and mean value. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 39: LST (above) and their histograms (below) for the proxy of solar radiance exposure (left) and 1 km LST (right) within the case study area. Slika 39: Temperature površja (zgoraj) in njuni histogrami (spodaj) za približek osončenosti (levo) in temperaturo površja v ločljivosti 1 km (desno) na testnem območju. 4.3.3. Data fusion by wavelet transform The Daubechies wavelet family (Chen and Zhu, 2007) of the fifth order was the optimal choice for the case study wavelet transform. Tests were also performed with other wavelet families however, the results were not as good. After both (low and high spatial resolution) datasets were transformed, the wavelet coefficients had to be appropriately fused together. The decomposition was performed up to level four (1 km LST). Only the 1 km LST coefficients were used for the reconstruction of the level four decomposition and all higher decomposition level coefficients were averaged by the weighted sum of corresponding wavelet coefficients: 0.7 and 0.3 for the third level (500 m), 0.5 and 0.5 for the second level (250 m) and 0.3 and 0.7 for the first level (125 m); the first weight stands for the 1 km LST while the second stands for the solar radiance exposure proxy. These weights yielded into the highest correlation with LST reference. The visual inspection of the downscaled (fused) image shows that this procedure works (figure 39) – one can observe that the low resolution data is appropriately fused with the high resolution data – for example, high temperatures in the left low corner (a plane that is not a part of the mountains) show the effect of low resolution data (figure 39 above right); on the other hand, the effect of the relief morphology Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. is visible all across the image and this is the consequence of the high resolution data (figure 40 left). Figure 40: Downscaled (125 m) LST (left) and the reference LST (right) in the case study area. Slika 40: Referenčna (desno) in temperatura površja z izboljšano ločljivostjo (125 m; levo) na testnem območju. The downscaled LST (figure 40 left) was compared to the Landsat TM 5 LST (figure 40 right). The extreme differences amounted to approximately ±7 °C with RMSE of 1.5 °C and the correlation index of 0.71. The results are quite similar to the statistics computed for the procedure of downscaling from 5000 m to 1000 m as described in the previous chapter. However, the visual comparison between the downscaled and reference data shows that the downscaled LST is overestimated regarding Landsat LST for the entire area, although a number of details in the downscaled image correspond to the reference image. There are many reasons that might yield in this systematic error. First of all, the reference LST was computed with the assumption that the Landsat TM 5 looks upon a black body – no emissivity data was used for reference LST computation, thus it might not be appropriately estimated. Secondly, linear scaling is a weak method of histogram matching. Certain errors might also occur due to input data – LST downscaled to 1000 m already has RMSE of approximately 2 °C, but the solar radiance exposure LST proxy should be estimated accurately enough. However, numerous details can also be seen in the flat terrain in the reference LST image (figure 40 right), which actually corresponds to a different land cover, and this was not taken into consideration. Figure 41 shows the differences between the downscaled and reference LST. Once again it can be seen that the downscaled LST is overestimated (the mean difference is approximately 1.5 °C). The differences are also great on the flat terrain (in the corner left Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. below), where the vegetation has a much greater influence on the LST distribution than the solar radiance exposure. Figure 41: Differences between the downscaled and reference LST. Slika 41: Razlike med temperaturo površja izboljšane ločljivosti in referenčnimi podatki. 4.4 Discussion Wavelet transform is a fast and easy to implement method for LST downscaling. It provides greater details in the results than, for example, contextual analyses, but the results of downscaling LST are not of the same quality as the results of panchromatic sharpening. The publications in the downscaling field are rare, only one similar study (regarding the used Landsat TM 5 data) was found (Venkateshwarlu et al., 2004), but the approach was tested in a case study area with completely different characteristics, thus it cannot be compared with the presented approach. The main problem of the presented approach is the histogram matching of the solar radiance exposure LST proxy data. The matching was based mathematically and not physically – the histogram was scaled and translated, but its shape was still not similar to the shape of the LST histogram. LST should be fully simulated in order to obtain greater similarity between both histograms. For example, in order to develop a new LST model and use it within the fusion procedure different geological, soil and vegetation parameters should be used alongside the relief morphology parameters. Even if the case study had not contained any flat relief (where the differences are the greatest), the histogram differences would have remained and the results would not have improved. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. One can conclude that downscaling LST from 1000 m to 125 m with only solar radiance exposure proxy using wavelets does not provide the expected results. The accuracy of approximately 1 °C was expected but not achieved (MODIS LST product has the accuracy of 1 °C already at the spatial resolution of 1000 m). However, if one considers that LST is a spatially inhomogeneous variable (for example, on a summer afternoon a road can be as much as ten degrees hotter than the grass merely a few meters away) the desired accuracy seems impossible to reach unless LST was modelled in a spatial resolution below 5 m (width of a road). Therefore, it might be concluded that solar radiance exposure has a great influence on LST, but the effect of the land cover or soil type can be even greater (at least in the areas that are not rough). This alternative hypothesis is also interesting at night, since the Sun is down and the solar radiance exposure is constant (zero) all across the area, yet LST varies across the area. Considering the unreliable results of the spatial resolution of 125 m, the 1000 m LST spatial resolution is a more appropriate input for AAT determination than LST. This is described in the following chapter: AAT is first determined from a 5000 m resolution LST, then from the 1000 m LST downscaled only with NDVI and on the end, AAT is determined also from the LST model described in chapter three. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. AMBIENT AIR TEMPERATURE PARAMETERIZATION FROM LAND SURFACE TEMPERATURE Ambient air temperature (AAT) can be interpolated with high accuracy into a continuous AAT field when sufficient measurements are available (chapter two). On the other hand, appropriate data for a continuous AAT field derivation (especially over large areas) can be obtained by remote sensing. However, the thermal signal of the air just above the surface is too weak when the signals of the surface or a higher atmosphere layer are taken into account, thus AAT can be measured directly by remote sensing from satellites only with an accurate AAT first guess. It is possible to observe the temperature at higher atmospheric layers by using a weighting function, which describes the change of the total signal transmittance with respect to the air-pressure that changes significantly at higher altitudes. On the other hand, the differences within the air pressure just a few meters above the ground are too small to reach an accurate AAT estimate (the ratio between the noise and the signal is then too high). Nevertheless, using additional information from the digital elevation model to model AAT from satellite sounder measurements was tested, but only in a low spatial (5000 m) and temporal (once per day) resolution (Méndez, 2004). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. AAT is not driven directly by the Sun but indirectly by the land surface temperature (LST). Therefore, the idea of using satellite images for parameterizing AAT from LST could be realized once the first LST data became available approximately 30 years ago. Numerous satellite sensors are appropriate for LST detection; geostationary satellites enable LST detection multiple times per day, for example the SEVIRI sensor on board MSG has a temporal resolution of 15 min and since most meteorological stations do not offer temporal resolution better than one hour, the use of this data provides an alternative to the available temperature diagnostic models. The temperature forecast in Slovenia is based on the ALADIN/SI model (ARSO, 2007b). The accuracy of the ALADIN/SI, ECMWF and other numerical weather prediction models and their correction is between 2 and 3 °C for a 24 h forecast (accuracy is generally better for temperature extremes) but greater deviations are possible especially in certain special meteorological conditions (for example low clouds) when the accuracy can be as much as three times worse. Similar accuracy is also obtained in Germany (figure 42), where the LME model is used for forecasting the weather (DWD, 2007), thus the goal of this study was to set up a parameterization based on remote sensing data that would provide an AAT field of high spatial (1000 m) and temporal (30 min) resolution with the expected accuracy of 2 °C. RMSE [°C] 3.2 3 2.8 2.6 2.4 2.2 2 1.8 1.6 Whole year: RMSE = 2.29 °C Jan Feb Mar Apr May Jun Jul Aug Sep Okt Nov Dec Month Figure 42: RMSE of AAT noon forecast in Germany during the year 2002; the line with circles presents the error for 24 h AAT forecast (2.3 °C for the entire year; DWD, 2007). Slika 42: RMSE opoldanskih temperatur zraka, napovedanih v Nemčiji v letu 2002; linija s pikami predstavlja napako za 24 urno napoved temperature zraka (2,3 °C za vse leto; DWD, 2007). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 5.1 Published methods Numerous approaches to determine the AAT from LST were tested in the past, but the connection between the surface, vegetation and air is so complex that the exact solution is still not available. The already tested approaches can be divided into four groups: simple statistical approach, advanced statistical approach, temperature-vegetation index (TVX) approach and energy-balance approach. Most of the described methods require LST measured by satellites, which is available only in sunny weather. A reader should also be aware that most of them function satisfactory only in certain conditions (for example no wind). The simple statistical approaches are usually based on the linear regression between LST and AAT. In certain conditions (no wind), they work excellently, but their quality is too poor for a general purpose. Chen et al. (1983) used the LST derived from a Visible and Infrared Spin Scan Radiometer (VISSR; 10.5–12.6 µm) data and compared it to the AAT on the height of 1.5 m. The tests were preformed in Florida, on cloudless nights, during four winters. They yielded a mean correlation index of 0.87 and an average sample standard deviation from regression of 1.6 °C. Davis and Tarpley (1983) estimated the daily minimum and minimum shelter temperatures for agricultural monitoring. For the NOAA-6 and NOAA-7 cloudless and partly cloudy retrievals the regression estimates based solely on temperature predictors from the satellite yielded residual standard deviations between 1.6 and 2.6 °C. Greena and Hay (2002) used AVHRR data for providing climatic variables across Africa and Europe for epidemiological applications. The correlations were examined between monthly mean AAT and variables derived from the AVHRR LST in Europe and Africa. Over a period of 12 months the mean accuracy of the retrieved AAT in Europe equalled 2.4 °C. The advanced statistical approaches consider at least two parameters. Cresswell (1999) estimated AAT from METOSAT LST by assuming that the solar zenith angle can be used as a proxy for the solar energy that reaches the surface. The AAT was then computed from the LST corrected by the function of the solar zenith angle. For approximately 72% of the processed temperatures an accuracy of 3 °C was achieved for all weather conditions in Angola, Namibia, South Africa, Zimbabwe and Botswana. EARS (2003) developed a methodology for the AAT estimation based on the comparison of noon and midnight Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. temperatures (in the case of a perfect heat exchange, there is no temperature difference between the surface and the atmosphere at night). The correlation index across the entire year for China equals 0.71. Jang (2004) used multi-layer feed-forward neural networks from the five AVHRR bands, elevation, solar zenith angle, and day of the year in order to estimate the AAT in southern Québec. The best accuracy was obtained with 22 hidden nodes that yielded to the correlation index of 0.93 and the RMSE of 1.8 °C. The temperature-vegetation index (TVX) approach is also a statistical approach but is described separately because it was tested and documented by numerous authors (e.g. Saravanapavan and Dye, 1995; Prihodko and Goward, 1997; Czajkowski et al., 2004). The approach is based on the hypothesis that the quantity of vegetation is correlated to LST (figure 43; less vegetation yields higher LST during the daytime). Furthermore, one can assume that the temperature on top of the canopy is the same as within the canopy if the canopy is infinitely thick – or to put it simply: in the event of a thick canopy LST equals AAT. The correlation between vegetation (NDVI) and LST is linear; the regression parameters are determined within the contextual (moving window) analysis. The approach has been mainly tested together with LST obtained from AVHRR data. The moving window size is usually 9 by 9 cells large. The approach was tested in different parts of the world and in most cases the RMSE is 3 to 4 °C for all weather conditions. Figure 43: TVX method: the regression line (right) parameters are determined within a moving window (light grey) for its centre cell. Slika 43: Metoda TVX: parametri regresijske premice so določeni za srednjo celico znotraj premikajočega se okna (svetlo modro). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The surface energy balance parameterizations are based on physics: the incoming net radiation flux with the anthropogenic heat flux equals the sum of the soil heat flux, the sensible flux and the latent heat flux (figure 44). An accurate modelling of these variables requires numerous data, thus a number of different approximations can yield a number of different models. Pape and Löffler (2004) performed a case study in the Norwegian mountains. The input data for the parameterization is mainly represented by ground measurements of different environmental parameters, but the authors expect that their parameterization can also be easily used with remote sensing data. The other required inputs are relief, substrate, vegetation type and calculations of the spatial distribution of the actual and maximum possible radiation depending on the relief. The deviations range between 0.2 and 2.1 °C and the correlation index equals 0.9. Sun et al. (2005) used the crop water stress index as the basis for establishing the relationship between air temperature and land surface temperature. Determining the crop water stress index and the aerodynamic resistance are crucial for this method. The crop water stress index can be determined from meteorological data or microwave remote sensing data. The aerodynamic roughness was determined according to land use type using an experimental look-up table. This method was applied for two MODIS datasets in a region located in the North China plain. An accuracy of 3 °C was achieved for over 80% of the processed data. There is also an operational parameterization within the Meteonorm software (Meteotest, 2004). Since the numerical coefficients were adjusted to the measured air temperatures on five Swiss meteorological stations it is a simplified surface energy balance parameterization. The parameterization considers: net radiation, cloud cover, wind speed and albedo but no water vapour or humidity are taken into consideration. For the daytime the correlation index equals approximately 0.87 (maximum deviation of 4.7 °C). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 44: Energy balance on the land surface. Slika 44: Energijsko ravnotežje površja. The purpose of this study was to choose an optimal approach that is accurate and simple for near real time operational use. Therefore, the methods that operate merely in certain special conditions were not considered. Furthermore, it is preferred if the parameterization is based on the boundary layer physics for this enable plenty possibilities for adjustments. Considering these two limitations, TVX and the surface energy balance methods were chosen for a detailed examination (table 7). It is clear that the TVX method is easy to implement, however it might provide poor results in the morning when the slope of the regression line can be close to zero or even positive (although it should be always negative according to the theory – a greater NDVI means a lower LST; the hypothesis cannot be confirmed for the mornings because overnight the bare ground cools down more than the vegetated ground). The Pape and Sun methods are based on the balance of energy between net radiation, soil heat flux, latent heat flux and sensible heat flux. They seem to have fewer limitations however they also require plenty of different inputs, which are often not available in near real time operational applications. The last method, described by Meteotest, seems to be simple and accurate. Furthermore, it has already been operational for some years and it has enough space for improvements, which makes it the optimal choice. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Table 7: Comparison of the characteristics between the potential methods to be used for AAT determination. Preglednica 7: Primerjava lastnosti najprimernejših tehnik za določitev temperature zraka. Author Input data Applicability Results quality Prihodko NVDI, LST AATSR, AVHRR, MODIS, SEVIRI r ~ 0.90 RMSE < 4.0 °C Pape incoming radiation, wind speed, relative humidity, precipitation, air pressure, elevation, soil type, etc. so far not tested with remote sensing data r > 0.90 . = 0.18–2.04 °C Sun LST, albedo, NVDI, landuse, air pressure, aerodyne. resistance, etc. AATSR, AVHRR, MODIS deviations: 0.18– 3.61 °C Meteotest LST, albedo, wind speed, incoming radiation so far not tested with remote sensing data r = 0.87 deviations < 4.7 °C 5.2 Testing existing methods The TVX (Prihodko and Goward, 1997) and Meteotest (2004) methods were tested in two case study areas: Slovenia and Germany. The characteristics of the Slovenian test were already described in chapter two (2.2). The relief in Germany is not as rough, except in the Alpine area in the south. Its position in Central Europe is limited with the Baltic and the North Sea in the North. The climate in Germany is influenced by a number of parameters but it changes gradually from the oceanic to continental and mountainous climate. 5.2.1. Data Most of the data used in this study was recorded by the Spinning Enhanced Visible and Infra- Red Imager (SEVIRI on board MSG) and processed by LandSAF (2007). At the moment, the used LandSAF products (LST, DSSF, DSLF, AL) are still in their pre-operational phase, which means that they should satisfy the major application requirements with certain limitations (i.e., scientifically checked but not validated, use of non operational data source, etc.). In Central Europe, this data has a spatial resolution of approximately 5–6 km. LST, one of LandSAF products, was already described in chapter three (3.1). The second used LandSAF product was the down-welling surface short-wave radiation flux (DSSF) which refers to the radiative energy in the wavelength interval of 0.3–4.0 µm that reaches the Earth's surface (per time and surface unit). This depends on the solar zenith angle, cloud coverage, atmospheric absorption and surface albedo. Furthermore, DSSF is strongly anti-correlated Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. with the observable top-of-atmosphere reflectance, thus DSSF that reaches the ground is lower in the presence of clouds (or even when the aerosol ratio in the atmosphere is high). Therefore, the top-of-atmosphere albedo is first computed from the observed directional reflectance values; then it serves as the most important piece of input information for the simple physical parameterisation of the radiation transfer in the cloud-atmosphere-surface system. DSSF is derived every 30 minutes from the three SEVIRI channels: VIS 0.6 µm, NIR 0.8 µm and SWIR 1.6 µm. The information on cloud coverage is obtained from SAFNWC. The information on the atmospheric water vapour content comes from the ECMWF numerical weather prediction model. The expected DSSF accuracy is greater than 110 Wm-2 in the event of cloudy weather and greater than 40 Wm-2 in cloud-free conditions. Down-welling surface long-wave radiation flux (DSLF) is a consequence of radiative energy in the 4–100 µm spectra. It is a result of atmospheric absorption, emission and scattering within the entire atmospheric column, thus it depends on the vertical profiles of temperature and gaseous absorbers. It can be determined by the radiation that originates from the shallow layer close to the surface (80% in the layer 0–500 m above the ground). It is derived every 30 min, but plotting its values over time shows that the values change only every third hour, which is probably a consequence of an error within the product generation process because DSLF changes continuously and not once per three hours (figure 45). Otherwise, the DSLF quality depends on the accuracy of both – cloudy pixel detection (again by SAFNWC) and atmospheric column characterisation (temperature and humidity). The relative DSLF accuracy is greater than 10%. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 45: DSLF for some of the chosen German meteorological stations (marked with different symbols) on 30.10.2005; notice how the values remain constant over a three hour long period. The line shows how DSLF changes during the day. Slika 45: Navzdol usmerjeno dolgovalovno sevanje za izbrane meteorološke postaje v Nemčiji (označene z različnimi simboli) za 30. 10. 2005; vrednosti so konstante za tri ure dolga obdobja. Linija nakazuje, kako se parameter približno spreminja prek dneva. The land surface albedo (AL) quantifies the part of the energy that is absorbed and transformed into heat and latent fluxes. It is important for determining weather conditions at the atmospheric boundary layer. AL is generated once per day (an iterative scheme allows the composition of the information from five days). Three short-wave channels (VIS 0.6 µm, NIR 0.8 µm and SWIR 1.6 µm) are used in order to generate the directional-hemispherical albedo at local solar noon. Moreover, broadband albedo is derived for the visible, near-infrared and total short-wave wavelength ranges. The information on cloud cover is obtained from SAFNWC while the dynamic information on the atmospheric pressure and water vapour content comes from the ECMWF numerical weather prediction model. Undetected clouds and systematic errors in the aerosol optical thickness values are the most important error sources. The estimated relative accuracy of AL is 10%. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 46: Position of the German meteorological stations used within the study – only those marked with circles and an ID number were used (data source ESRI maps and GLOBE DEM – Hastings and Dunbar, 1998). Slika 46: Lega nemških meteoroloških postaj, uporabljenih v študiji – samo postaje, označene s krožci in identifikacijsko številko, so bile uporabljene na koncu (vir podatkov ESRI maps and GLOBE DEM – Hastings and Dunbar, 1998). For verification, the modelled AAT values were compared with the AAT measurements. Alongside the measurements provided by the Environmental Agency of the Republic of Slovenia, the meteorological measurements for Germany were also available (hourly were provided by the German Weather Service). The measurements from 49 stations were available, but some of them did not cover the entire year 2005, thus they were not included Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. into the database (33 stations). Some of these stations are situated too far north, where the SEVIRI data is no longer reliable, thus only 13 stations were used in the end (figure 46). 5.2.2. Testing Temperature Vegetation Index method Table 7 proposes to test four methods that use remote sensing data. It was decided to test only TVX (Saravanapavan and Dye, 1995; Prihodko and Goward, 1997; Czajkowski et al. 2004) and and Meteotest’s (2004) parameterizations, since the second and third parameterization require too much input data for near-real time applications. The testing of the existing methods was performed without the downscaling of the original data. Only data from cloudless or partially cloudy daytime weather was used – the list of the selected dates used for testing can be seen in Appendix D. The Temperature Vegetation Index (TVX) method was tested only in Slovenia; at first MODIS LST (1000 m) and NDVI (16-day composite, 1000 m spatial resolution) were used. Tests (for randomly selected partially cloudy or cloudless days – in the daytime) have shown that results of similar quality to the ones previously reported can be obtained: RMSE of 2.5 °C and r of 0.93 (figure 47) is even better than in most previous studies (Saravanapavan and Dye, 1995; Prihodko and Goward, 1997; Czajkowski et al. 2004). In order to improve the method three different moving window sizes were tested: 7 by 7, 9 by 9 and 11 by 11. The quality improved with the size of the moving window. The best results were obtained by using an 11 by 11 cells moving window with the NDVI of the “infinitely thick” canopy set to 0.86. The second test was performed by using the SEVIRI LST and a generalized MODIS NDVI. The results were similar to the ones obtained in the previous test (RMSE of 2.6 °C and r of 0.96; figure 48), which is again comparable to previous studies (the only problem might be a low spatial resolution). In this case the moving window was 7 by 7 cells; larger moving windows provide too generalized regression parameters (considering the 5 km spatial resolution), while smaller moving windows fail to provide sufficient data pairs to establish a reliable dependency between the two variables. All previous studies are based on AVHRR images working with a 9 by 9 moving window. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 47: Measured-parameterized AAT scatterplot; MODIS LST and NDVI; r = 0.93; RMSE = 2.5 °C. Slika 47: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka; temperatura površja, pridobljena in NDVI z MODIS-om; r = 0,93; RMSE = 2,5 °C. Figure 48: Measured-parameterized AAT scatterplot; SEVIRI LST &MODIS NDVI; r = 0.96; RMSE = 2.6 °C. Slika 48: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka; temperatura površja, pridobljena s SEVIRI-jem, NDVI z MODIS-om; r = 0,96; RMSE = 2,6 °C. The MODIS vegetation index also consists of the enhanced vegetation index (EVI), thus the method was also tested by using EVI instead of NDVI. Within areas densely covered by vegetation, EVI has an advantage over NDVI – NDVI is quickly saturated and provides no “contrast” in densely vegetated areas. On the other hand, while EVI provides a greater contrast in the densely vegetated areas, its performance is less satisfactory in areas with low Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. vegetation density. Since more than 60% of Slovenia is covered with forest, EVI was tested as a proxy for vegetation instead of NDVI. The results were similar in highly vegetated areas but worse than the NDVI results in urban areas – the correlation between EVI and LST was weaker than the correlation between NDVI and LST (figure 49). Figure 49: Relationship NDVI-LST and EVI-LST on the 27.10.2005 south of Ljubljana (Slovenia; Zidarič, 2007). Slika 49: Povezava med NDVI ter EVI s temperaturo površja za 27. 10. 2005 južno od Ljubljane (Zidarič, 2007). However, the TVX method woks satisfactory only during the daytime – the slopes computed by contextual analyses are often positive, which is impossible according to the background theory (with a larger NDVI the LST should be smaller and not vice versa). Positive slopes in the regression line are more often found in the morning (this can be observed if one uses the SEVIRI LST, which is available every 15 min) than in the afternoon. This occurs because the surface cools down during the night – the bare ground usually cools down faster and to a greater extent than a vegetated surface, thus the TVX hypothesis is not valid in the morning. The positive slopes of the regression line also have a mathematical background. Some areas are covered in vegetation and the NDVI is high (more than 0.8) with only slight variability (homogeny). In these areas the LST is usually constant, thus the results of regression analysis can change significantly even with a small change in the data – the results are statistically uncertain, because the temperature variability in these cases cannot be explained merely with NDVI. Therefore, one can conclude that it makes greater sense to equalize the LST values with AAT in the homogeneous areas with a high NDVI. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The method seems easy to implement, however it has certain limitations. First of all, it is almost impossible to obtain an AAT of high temporal and spatial resolution. By using the SEVIRI LST, one can expect a high temporal resolution; however the spatial resolution can be relatively high only in Africa, where the sensor records LST in nadir (spatial resolution of 3.1 km). It is theoretically possible to have a high temporal and spatial resolution close to the Earth poles, for a number of satellites have a near polar orbit and cross the pole areas once every one to two hours. The second limitation of the TVX method is the problem with the positive slope of the regression line. Due to the snow, which has a low NDVI, the winters might cause an even greater problem and place the entire TVX hypothesis under question. 5.2.3. Testing Meteotest method The parameterization used by Meteotest seemed especially interesting because it is already implemented in their software (Meteonorm). The numerical coefficients of the parameterization were fitted to merely five Swiss meteorological stations (correlation index of 0.87 for daytime, with the maximum deviation of 4.7 °C), thus the parameterization might not provide satisfactory accuracy in other areas. Two different parameterizations were derived for AAT. The daytime parameterization is used when DSSF is greater than 5 Wm-2. The parameterization also considers the wind velocity u. AAT =LST -(0.0015·(1- AL)· DSSF - 0.7)·exp(- 0.09u) (22) Additional conditions are required during the night time (if AAT is less than -3 °C, then the results are valid only when LST is less than 0 °C; N presents cloud coverage in oktas, where one okta corresponds to one eighth of the sky). 32 AAT =LST -(0.0006· N - 0.037 · N + 0.376 · N - 4.7)·exp(- 0.218u) (23) The parameterization was first tested in Slovenia for the period between 30.05.2005 and 24.12.2005. Only daytime data from cloudless or partially cloudy days was used – the list of the selected dates used for testing can be seen in Appendix D. The first test showed that there is a high correlation between the measured temperatures and the temperatures parameterized with the Meteotest equation, although the results did not fit perfectly (figure 50). In order to improve the parameterization, the correlation between AAT and the remote sensing data and Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. meteorological variables was estimated. As expected the strongest correlation was found between AAT and LST. AAT is also strongly correlated with DSSF. The differences between LST and AAT are the highest at the daily peak of DSSF. Furthermore, these differences are smaller in the presence of a strong wind. Figure 50: Measured-modelled AAT scatterplot for selected days at 6:00, 13:00 and 20:00 GMT; results obtained in Slovenia with the original Meteotest parameterization; r = 0.91; RMSE = 4.0 °C. Slika 50: Razsevni grafikon med merjeno in modelirano temperaturo zraka za izbrane dni ob 6:00, 13:00 in 20:00 GMT za izbrane postaje v Sloveniji; r = 0,91; RMSE = 4,0 °C. The Meteotest parameterization for the night time is based on cloud coverage derived from LandSAF data. This was computed from the no-data values in the LST layer within a 3 by 3 window. The spatial resolution of the SEVIRI is already relatively low in Central Europe and the window of this size covers an area of 15 by 15 km, which is large enough to estimate cloud coverage. On the other hand, the estimation is made from nine cells only and usually all of them are covered with clouds or all of them are cloudless, thus the estimated cloud coverage has binary characteristics, which led to an unsatisfactory accuracy. Since cloud coverage can be understood as a proxy for the radiation coming from the atmosphere, LandSAF DSLF was tested as the input data instead of cloud coverage in oktas. The results were encouraging and thus another parameterization was generated, this time one that works for night and daytime. The coefficients were calibrated by the measurements from 12 stations in Slovenia. AAT= 1.13·(LST -(0.012 ·(1- AL)· DSSF - 0.008· DSLF )·exp(- 0.29u))+1.65 (24) Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The results were much better as regards the original parameterization. However, there were still some residuals larger than ±7 °C. It was assumed that the low deviations were a consequence of clouds (sub-pixel size). Furthermore, there is a high probability that these sub- pixel clouds exist next to the pixels that are already classified as clouds. Therefore, if one of the four closest pixels (AAT is estimated through a bilinear interpolation from the four closest pixels for each dataset and if one of them is “contaminated”, the interpolated value is an error source) to the selected meteorological station was covered with clouds, then the AAT was not computed for that station. The same parameterization as above was used once more and it seemed to be robust because the best fit line did not change significantly, yet the quality of the parameterized AAT improved (figure 51). Figure 51: Measured-modelled AAT scatterplot for selected days at 6:00, 13:00 and 20:00 GMT; results obtained in Slovenia with a changed Meteotest parameterization; r = 0.97; RMSE = 2.3 °C. Slika 51: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane ob 6:00, 13:00 in 20:00 GMT za izbrane postaje v Sloveniji; r = 0,97; RMSE = 2,3 °C. The residuals were also validated for every station. In general the residuals for station 38 (situated in the Alps; figure 3) are much higher than zero (the measured AAT is lower than the estimated AAT). The reason for this might be found in the relief of the area – the station is situated on a slope facing north-west, while its surrounding area is situated on a huge slope facing south. The area with the station does not receive as much solar irradiation as the surrounding area, thus the station can not be representative for the surrounding area. The problems of representativity also occurred in the study area of Germany. The first tests Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. showed that using a Meteotest based parameterization defined for Slovenia does not function as accurately in Germany as it does in Slovenia. The scatterplot (figure 52) shows residuals for certain measurements at 20:00; the marked points within the ellipse are measurements made by the Zugspitze meteorological station, which is obviously not representative for the surrounding area (because of its elevation), thus it was removed from the further analysis and the result improved slightly (figure 53 to figure 56). Figure 52: Measured-modelled AAT scatterplot for selected days at 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.75; RMSE = 4.5 °C; the points within the ellipse belong to the Zugspitze meteorological station. Slika 52: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane ob 20:00 GMT za izbrane postaje v Nemčiji; r = 0,75; RMSE = 4,5 °C; točke znotraj elipse pripadajo meteorološki postaji Zugspitze. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 53: Measured-modelled AAT scatterplot for selected days at 6:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.90; RMSE = 3.3 °C. Slika 53: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane ob 6:00 GMT za izbrane postaje v Nemčiji; r = 0,90; RMSE = 3,3 °C. Figure 54: Measured-modelled AAT scatterplot for selected days at 13:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.88; RMSE = 6.4 °C. Slika 54: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane ob 13:00 GMT za izbrane postaje v Nemčiji; r = 0,88; RMSE = 6,4 °C. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 55: Measured-modelled AAT scatterplot for selected days at 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.96; RMSE = 2.9 °C. Slika 55: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane ob 20:00 GMT za izbrane postaje v Nemčiji; r = 0,96; RMSE = 2,9 °C. Figure 56: Measured-modelled AAT scatterplot for selected days at 6:00, 13:00 or 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.89; RMSE = 4.4 °C. Slika 56: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane ob 6:00, 13:00 ali 20:00 GMT za izbrane postaje v Nemčiji; r = 0,89; RMSE = 4,4 °C. The estimated accuracy of the thus modelled AAT was still far below the desired accuracy (RMSE of 4.4 °C and correlation index of 0.89 is meagre when taking the desired 2 °C into account). There was no significant change even when the coefficients from formula (24) were changed for a better fit. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. AAT= 1.13(LST -(0.012 ·(1- AL)· DSSF - 0.008· DSLF )·exp(- 0.12u))+1.00 (25) Additional LandSAF data was used because the measurements for Germany were available for every hour, which made it possible to observe how the residuals change during the day. The scatterplot (figure 57) shows the distribution of all available data using the equation (25). Figure 57: Measured-modelled AAT scatterplot for selected days from 4:00 to 20:00 GMT; results obtained in Germany with a changed Meteotest parameterization; r = 0.90; RMSE = 3.8 °C. Slika 57: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane od 4:00 do 20:00 GMT za izbrane postaje v Nemčiji r = 0,90; RMSE = 3,8 °C. 5.3 A new AAT parameterization The results using the Meteotest based parameterization in Slovenia are more promising than in Germany, however there are certain limitations in combining the used parameterization and the used data, thus the higher quality of results was reached only after certain changes were applied into the procedure. The scatterplots (figure 53 to figure 55) from different times of the day suggest that there might be a daily periodical parameter that influences the used parameterization of LST into AAT. Therefore, the data that ranges between 4:00 and 20:00 GMT on the dates between 03.07.2005 and 12.11.2005 (previous scatterplot – figure 57) was used to plot the distribution of the residuals (parameterized minus the measured AAT) as a function of time for every date with available data in order to recognize any significant pattern Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. that might help to improve the parameterization (figure 58). It was discovered that the residuals are generally positive during the daytime and the greatest around local noon. Within the analysis, only daytime data from cloudless or partially cloudy days was used – the list of the selected dates used for testing can be seen in Appendix D. Figure 58: Residuals (between the parameterized and measured AAT) [°C] as a function of daily time [h] (03.06.2005); each station has its own symbol. Slika 58: Odstopanja (med parameterizirano in merjeno temperaturo zraka) [°C] kot funkcija dnevnega časa[h] (03. 06. 2005) ; vsaka postaja je označena s svojim simbolom. 5.3.1. Differences between LST and AAT as a function of accumulated solar irradiation Figure 58 shows a pattern that is most likely a consequence of a turbulent heat flux, which can also be shown with the differences between LST and AAT. Therefore, another test was preformed in order to investigate these differences during the day. First a set of plots was performed with the residuals as a function of the accumulated solar irradiation, which is the sum of all solar energy received by the surface from sunrise until the chosen time. The accumulated solar irradiation has a strong influence on AAT and LST because both of these parameters not only depend on the current weather conditions but also on the weather in the past – it is warmer when more solar irradiation is received. As seen in the lines representing the residuals for each day (figure 59), the accumulated solar irradiation varies significantly Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. through the year; the Sun is for just a few hours in the sky during the winter and the accumulated solar irradiation is smaller than during the summer. Figure 59: Lines representing the residuals [°C] plotted as a function of accumulated solar irradiation [Whm-2] for selected days in 2005; more curved lines are made for winter days and less curved for summer days. Slika 59: Prilegajoče se linije odstopanj [°C], izrisana kot funkcija kumulativnega obseva [Whm-2] za izbrane dni v letu 2005; bolj ukrivljene črte so bile narejene na podlagi zimskih in manj ukrivljene na podlagi poletnih dni. 5.3.2. Differences between LST and AAT as a function of time Another set of plots was performed as a function of time, for this is easier to consider within the parameterization than the accumulated solar irradiation. As expected, the differences between LST and AAT are the greatest around noon and close to zero around sunrise and sunset (figure 60). During the summer the amplitude residuals are smaller and positive for a longer time when compared to the winter (figure 61). Furthermore, the plots have shown that the residuals are distributed more systematically over time than expected – the pattern in the plots of differences between LST and AAT is stronger than the pattern of residuals between the parameterized (changed Meteotest parameterization) and measured AAT. Therefore, it made no sense to correct the Meteotest based parameterization; instead it made greater sense to make a new parameterization as a function of time considering the daily and annual cycle. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 60: Differences between LST and AAT [°C] as a function of time of day [h] (03.06.2005); each station has its own symbol. Slika 60: Razlika med temperaturo površja in temperaturo zraka [°C] kot funkcija dnevnega časa [h] (03. 06. 2005); vsaka postaja je označena s svojim simbolom. Figure 61: Lines representing the residuals [°C] plotted as a function of the time of day [h] for a few days in 2005. Slika 61: Prilegajoče se linije odstopanj [°C], izrisana kot funkcija dnevnega časa [h] za izbrane dni v letu 2005. 5.3.3. Parameterization of the differences between LST and AAT The lines depicting the differences between LST and AAT (figure 61) have a shape of a sinus function reaching the amplitude around noon, which corresponds to the change in the zenith Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. angle of the sun. Therefore, the LST-AAT differences were first parameterized with the cosine of the zenith angle of the sun. Since the amplitude of this curve might be locally dependent, the cosine was multiplied with the natural logarithm of NDVI (in order to expand the range of NDVI, which was generally rather small in the case study area). The correlation between the real and parameterized differences reached almost 0.8, which is comparable to the similar approach used by Cresswell (1999). Furthermore, the Meteotest parameterization showed which parameters have a strong influence on the LST-AAT differences, thus LandSAF data was integrated into the parameterization. Multiple regression was used to calibrate the optimal coefficients (from 700 randomly selected measurements). Furthermore, LST had to be downscaled to a higher resolution (problems with mountainous stations in Slovenia and Germany), which made the adjustments to the local anomalies possible. A simple downscaling procedure was first applied to LST based on the correlation between NDVI and LST (the TVX approach hypothesis described in chapter 3.3; Figure 21). The final results are shown in figure 62, where all measurements (1473) are presented on the scatterplot. AAT=LST - 5.399 - 6.581·cos(z)·ln(NDVI )+ 0.032· DLSF - (26) - 0.014·(1- AL) · DSSF - 3.499·exp(- 0.3u) Figure 62: Measured-modelled AAT scatterplot for selected days between 4:00 and 20:00 GMT; results obtained in Germany with a new parameterization; r = 0.96; RMSE = 2.1 °C. Slika 62: Razsevni grafikon med merjeno in na novo parametrizirano temperaturo zraka za izbrane dni od 4:00 do 20:00 GMT za izbrane postaje v Nemčiji; r = 0, 96; RMSE = 2,1 °C. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The extreme residuals still remain quite high (minimum -8.0 °C, maximum 7.6 °C) but the overall results seem promising (RMSE = 2.1 °C, r of 0.96). The same parameterization was also used in Slovenia, but the results were slightly worse with RMSE reaching 2.7 °C. The results improved after the parameterization coefficients were changed as in (27): RMSE of 1.8 °C, r of 0.98, minimum of -6.6 °C and maximum residual of 5.5 °C. AAT=LST - 4.25 -1.27 ·cos(z)·ln(NDVI )+ 0.022· DLSF - - 0.0079 ·(1- AL) · DSSF - 2.99 ·exp(- 0.3u) (27) Figure 63 shows how RMSE varies over time (selected days are listed in appendix D). It can be observed that the worst accuracy is obtained at 6:00 GMT during the last three months in the year. The cause for this might be that some areas were covered by snow for a short while, because the albedo at some stations varies significantly over a short time. 0 0.5 1 1.5 2 2.5 3 3.5 4 RMSE [°C] 6:00 13:00 20:00 June July September Oct.–Dec. Figure 63: RMSE variation for selected situations. Slika 63: Variabilnost RMSE za izbrane situacije. As it can be observed by comparing (26) and (27), the coefficients by the cosine of the solar zenith angle differentiate significantly and this is the reason that the parameterization for one area is not accurate also in the other area. This parameter is usually highly correlated with DSSF, thus only one of them should be used to satisfy the mathematical rule about using independent data in multiple regression. However, the solar zenith angle can provide more information in the late afternoon than DSSF (after sunset when DSSF is already zero). On the Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. other hand, the solar zenith angle does not provide any information on the transparency in the atmosphere, which is contained in DSSF, thus it was decided to use both parameters. Wind velocity was the only parameterization input that was not provided as raster data within the study although this data exists in the spatial resolution of 2500 m and temporal resolution of 3 hours (ARSO, 2007a). However, wind is a parameter that has the least systematic effect on the presented parameterization. Therefore, another parameterization was performed for Slovenia without considering the wind velocity: RMSE of 1.9 °C, r of 0.97, minimum of -7.1 °C and maximum residual of 4.4 °C; figure 64 shows how RMSE varies over time (selected days are listed in appendix D). AAT=LST -6.634-1.434 ·cos() ( NDVI )+0.021·DLSF -0.0069 ·(1-AL z ·ln ) ·DSSF (28) 4 3.5 3 2.5 2 1.5 1 6:00 0.5 13:00 0 20:00 Figure 64: RMSE variation for selected situations (wind velocity not considered). Slika 64: Variabilnost RMSE za izbrane situacije (hitrost vetra ni bila upoštevana). This kind of parameterization offers the opportunity to finally compute a continuous AAT field (figure 65). LandSAF provides the required input data (DSSF and DSLF with a 30 min temporal resolution), thus the final results have a high temporal and spatial (1000 m due to the LST downscaling) resolution. RMSE [°C] June July September Oct.–Dec. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 65: An AAT example for Slovenia, generated with the described procedure (27.07.2005 at 20:00 GMT) in 1000 m spatial resolution. One can observe numerous no-data areas, which are the consequence of cloud coverage. Slika 65: Temperatura zraka, določena z opisanim postopkom za Slovenijo (primer za 27.07.2005 ob 20:00 GMT) v prostorski ločljivosti 1000 m. Videti je lahko tudi izbrane območij brez podatkov v belem, ki so posledica oblakov. Figure 66 shows a comparison between the residuals histogram obtained with or without considering wind velocity in the parameterization of LST into AAT. The comparison is important because it can help us understand the effect of the wind velocity on AAT parameterization. It can be seen that the residuals obtained in the parameterization considering wind velocity (27) are closer to zero than in the parameterization not considering wind velocity (28). Therefore, although the correlation between AAT and wind velocity is low, wind does influence AAT and it would be appropriate to include wind data into the parameterization. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 0.0 10.2 19.1 35.2 23.0 10.2 2.3 0.0 11.3 21.9 29.7 25.4 9.8 2.0 0 5 10 15 20 25 30 35 40 below -5 -5 to -2.5 -2.5 to -1 -1 to 1 1 to 2.5 2.5 to 5 abov e 5 Residuals [°C] Frequency [%] Figure 66: Histogram of residuals considering wind (blue) and not considering wind (red) in the parameterization. Slika 66: Histogram odstopanj v parametrizaciji z vetrom (modro) in brez njega (rdeče). 5.3.4. Parameterization of differences between LST and AAT form the LST downscaled with contextual analyses LST used in equations (26) to (28) was downscaled with a simple procedure based on the relationship between LST and NDVI. However, as described in chapter three further influences can be considered by downscaling. Therefore, the final step of this study was to examine if the parameterization can be additionally improved by LST containing more details. The procedure of downscaling LST was only tested in Slovenia, thus the AAT from the improved LST was determined for the same area. Surprisingly the results did not improve. Figure 67 shows numerous details due to the numerous details in the LST image, however the statistics got worse – the extreme residuals remained similar (minimum -6.4 °C and maximum 6.2 °C) while RMSE (2.2 °C) and the correlation index (0.97) got worse. Once more, the wind velocity was not used in the parameterization. AAT=LST +4.65-3.65·cos()z +0.013·DLSF -0.010 ·(1-AL)·DSSF (29) Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 67: AAT for Slovenia (an example for 28.07.2005 at 13:00 GMT) in 1000 m spatial resolution. Slika 67: Temperatura zraka za Slovenijo (primer za 28. 07. 2005 ob 13:00 GMT) v prostorski ločljivosti 1000 m. A possible reason for the poor statistics might lie in the fact that insufficient data was used to fit SEVIRI to MODIS LST – most MODIS data used for the fit was acquired during daytime by Aqua. Usually Aqua passes the equator around 13:00 GMT, thus the accuracy was tested also just for the AAT measurements at 13:00. The statistics improved – the correlation index and RMSE did not improve significantly comparing to the results acquired by the procedure, in which LST was downscaled only by NDVI, but the minimum of -3.7 °C and maximum residual of 5.4 °C (figure 68) represent a significant improvement. Therefore, the LST downscaled as presented in chapter three should be used in the morning or at midday, but otherwise it is better to use the simplified downscaling procedure based on NDVI without SEVIRI LST correction. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 68: Measured-modelled AAT scatterplot for selected days at 13:00 GMT; results obtained in Slovenia with the use of the improved downscaling procedure on a few stations; r = 0.97; RMSE = 1.8 °C. Slika 68: Razsevni grafikon med merjeno in parametrizirano temperaturo zraka za izbrane ob 13:00 GMT za izbrane postaje v Sloveniji (uporabljen izboljšan postopek izboljšave ločljivosti temperature površja); r = 0, 97; RMSE = 1.8 °C. 5.4 Discussion The presented approach makes it possible to estimate a continuous AAT field of high temporal and spatial resolution from the data that is easily obtained from the internet even free of charge. Therefore, it provides an excellent alternative to the previous approaches that might require certain input data that is difficult to obtain. The LandSAF products are still not in the fully operational phase, some product validation processes are finished but there might still be some unsolved problems as shown in the DSLF description (this parameter changes gradually over the day, thus it does not have such a large affect on the overall quality of the presented parameterization). The MODIS products have already been operational for a few years and one can expect them to be more reliable. MODIS LST was used only to test the TVX approach and MODIS NDVI was additionally introduced into the parameterization for the downscaling procedure. Additional data that will need greater attention in the future is wind data. Only wind velocity measurements at meteorological stations were used within this study, thus the aim of the Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. study cannot be fulfilled – a continuous wind velocity input is needed; the wind analysis field from ECMWF for example can be implemented in the future. ECMWF provides wind velocity at the height of 10 m and a spatial resolution of 0.1 angular degree every four hours and even better data is available for Slovenia – DADA provides wind velocity data in the spatial resolution of 2500 m and temporal resolution of 3 hours (ARSO, 2007a). This data should be adequate for the proposed parameterization. The downscaling procedure used in the study improves the parameterization quality, which means that spatial resolution is important. There are numerous meteorological circumstances with a strong local effect (thermal inversion, rough relief, heat islands in cities, etc.) that can be presented only with high spatial resolution data. During the daytime the quality of the AAT obtained with the presented parameterization on the selected cloudless or partially cloudy days is better than the average quality of a 24 h AAT forecast (ARSO, 2007b; DWD, 2007). The results should be actually compared with diagnostic models, but the information about their quality is usually unknown. Furthermore, the quality of the presented parameterization is comparable or even better than the quality of most previous approaches that use remote sensing data for AAT parameterization: the results are similar to the results of Sun et al. (2005; most recent paper in this field) and much better than any results obtained with the TVX approach. However, the results are a little worse than Jang (2004), but it has to be understood that all these approaches are based on AVHRR or MODIS data, which have a higher spatial but poorer temporal resolution. The drawback of the presented parameterization is its regional dependency – the numerical coefficients differentiate for Slovenia and Germany. The reason for these differences is still unknown, thus it should be investigated in the future; the reason might be revealed when the parameterization is tested in at least one more case study area. One of the possible reasons for the differences can be satellite image distortions that occur due to the MSG-Earth geometry – the view angle of SEVIRI (on board MSG) changes significantly within Europe, which yields into inconstant spatial resolution. The reason might also be that some other parameter was not taken into consideration, thus it makes sense to check the parameterization in a third independent area. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. CONCLUSION During his lifetime my fellow countryman Herman Potočnik “Noordung” was considered to be a dreamer because of his descriptions of space travel that he described in his book “Das Problem der Befahrung des Weltraums”. However, a lot of his visions came true, even geostationary satellites, the orbit of which he presented in 1928. I guess he never thought that, eighty years later, these satellites will be used to measure the land surface temperature (LST) and indirectly determine the ambient air temperature (AAT). Even though my thesis cannot be compared to the life-work of Potočnik, it presents some interesting results: • First the SEVIRI LST was downscaled with the contextual analysis from 5000 to 1000 m and then to 125 m with wavelet transform (the LST model in 1000 m is comparable to the MODIS LST, which was considered a reference LST, but the quality of the 125 m LST is too poor to be implemented into the further procedure – AAT determination). • AAT was computed by using two approaches – by interpolating ground measurements and by modelling from LST (when AAT measurements are available, the interpolation is optimal solution since it provides better results). For the conclusion a comparison between both approaches of determining the continuous AAT field was performed. Both approaches are based on multiple regression, thus the computational speed is similar in both cases. The mean absolute difference of the results obtained from the meteorological stations equals 1.2 °C with the standard deviation of 1.5 °C. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 69: Bias (above) and the accuracy estimation for both approaches of AAT determination. Slika 69: Sistematično odstopanje (zgoraj) in ocena točnosti obeh pristopov izračuna temperature zraka. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Figure 69 shows the spatial distribution of the accuracy parameters for both approaches for selected days at the selected meteorological stations. It can be seen that small biases (mean residual) correspond well to high accuracy and vice versa. The accuracy obtained by the interpolation is generally better than the accuracy obtained by the parameterization from LST. In general, the accuracy is better in the eastern part of Slovenia, where the relief is flat. The wind velocity data was not considered in the interpolation approach, nor in the final version of the parameterization from LST, thus a sensitivity analysis was preformed in order to estimate the influence the wind velocity has on the accuracy of the results. Situations with wind velocity below 1 m/s were considered as calm and situations with wind velocity greater than 4 m/s were considered as windy. The interpolation approach seems very robust because the absolute residual values are uncorrelated with the wind velocity (r = -0.04). Moreover, the accuracy in windy situations is even slightly better than in calm situations. The reason for such robustness is probably the use of relief derivates such as relief slope or TduP, which have a good correlation with the air masses movement. On the other hand, the accuracy of AAT parameterized from LST in windy situations is approximately 25% worse than the accuracy in calm situations, but the correlation is still low (r = -0.14). The wind velocity of over 20 m/s is not unusual in western Slovenia (it was always below 5 m/s within the selected cases), thus one can assume that the accuracy in such cases is even worse. Therefore, AAT parameterized from LST can be successfully determined only in relatively calm situations (in the end, it was determined that the accuracy is insensitive to wind velocity below 3 m/s while stronger wind has a strong influence on the results). Another limitation of the AAT parameterization from LST is snow, because the ground albedo rises and the parameterized temperatures are underestimated. To sum up, four hypotheses were examined; the first one was rejected – i.e. modelling LST in a high spatial and temporal resolution with a similar accuracy as the reference data. The accuracy of the downscaled LST to 1000 m is twice as bad as the accuracy of the reference MODIS LST, which has the accuracy of 1 °C. Moreover, the relationship between the solar radiance exposure and LST exists, but it is too weak to use in the process of downscaling LST to 125 m spatial resolution in areas with inhomogeneous land cover. Therefore, it is not possible to model LST at 125 m spatial resolution but it is possible to model LST at 1000 m Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. spatial resolution, with 15 min temporal resolution, and the accuracy of approximately 2 °C, which is comparable to the original SEVIRI LST accuracy. The second hypothesis was accepted – i.e. modelling the continuous AAT field from a high spatial and temporal resolution LST model with a higher accuracy than the accuracy of the meteorological forecast. The proposed model makes it possible to model the AAT for a cloudless sky during daytime with the accuracy of 2.1 °C, which is slightly better than the 24 h meteorological forecast accuracy for all weather situations. Moreover, the proposed procedure achieves a similar or even better accuracy than the previous similar studies with a better temporal and spatial accuracy. The third hypothesis was accepted – i.e. the accuracy of the AAT parameterized from remote sensing data is better when the weather is calm and still. The accuracy in the event that the wind velocity is higher than 4 m/s is 25% worse than the accuracy in calm situations. Therefore, the AAT parameterization is appropriate only in calm situations (wind can often be stronger, e.g. in the western part of Slovenia). The fourth hypothesis was accepted – i.e. using stationary data for the interpolation of AAT into a two-dimensional field. It provided results with the accuracy of 1.6 °C, which is superior to the accuracy of AAT computed from LST. The reason is probably higher spatial resolution, which can also take into consideration the local anomalies. However, the interpolation is based on AAT measurements, which are not always available. The topics investigated in the thesis have not been intensively researched in the past. Certain examples of using LST to compute AAT can be found in literature, but most of them deal merely with special situations because the problem seems to be too complex for a universal solution. The situation is even more difficult in the downscaling field – wavelets transform, which has been often described during the past years seems a good solution for fusing similar data, but it is not a good method for LST downscaling. Only one article dealing with downscaling LST was found, which means that there remains a lot to research. To conclude, numerous different methods for statistical modelling of AAT from a number of predictors were tested in the thesis. This might appear confusing sometimes, because the connection between some chapters might not be clear at first sight. The purpose of the thesis Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. was to show that one needs a continuous AAT field and if there are sufficient measurements they should be interpolated. In the contrary case, using remote sensing data is an appropriate solution but the results yield poorer accuracy compared to interpolation. It also has to be understood that both approaches might produce some no-data areas – the first method in areas without appropriate data and the second one in areas covered by clouds. Due to this disadvantage as regards the clouds, it is most likely that the remote sensing approach represents the most interesting one in the solar energy field, where AAT is needed especially when the Sun shines on the solar devices. A “side effect” of the described approach is an LST model of high spatial and temporal resolution that was tested in Slovenia. It can be used in numerous applications regarding the interactions between soil, vegetation and the atmosphere. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. SUMMARY The goal of the thesis is to show the methodology that can generate a continuous ambient air temperature (AAT) field of a high spatial and temporal resolution. The interpolation of spatial data recorded in meteorological stations is difficult, especially in mountainous countries with low autocorrelation among measurements. As an interpolation method that minimizes the interpolation error kriging is not adapted to inhomogeneous areas, thus a statistical method based on multiple regression was used. Temperature measurements on a set of 20 synoptic stations in Slovenia were made three times a day (7:00, 14:00 and 21:00 Central European Time) throughout 2005, which means that there are 1095 situations to be interpolated in total. The measurements were linked with the data stored in GIS: firstly, with the DEM with derived layers such as relief slope, relief aspect, etc., and secondly with the land cover including derived layers such as distance from forest. Following the first test, the results were not satisfactory due to numerous high residuals. Step by step, different tests were performed in order to improve the results: first, the distance from the sea was added as an explanatory variable; second, the distance from the forest was removed. Finally, MODIS images were introduced into the GIS and used within the multiple regression. The results improved after each test, with the standard deviation progressively improving to the final 1.6 °C. However, sometimes the AAT measurements are too sparse (spatially or temporally) to determine an AAT field of sufficient quality for the entire area of interest (at a given moment) especially in the case of inhomogeneous areas. Some applications require AAT data in high spatio-temporal resolution, partly even in near-real-time. Thus a parameterization, which determines the AAT from remote sensing data, was also developed. Since the ambient air Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. temperature is not driven directly by the Sun, but indirectly by the ground, a number of previous studies used land surface temperature (LST) to parameterize AAT. The present study is based on the multiple regression of the SEVIRI data (on board MSG). The input parameters – land surface temperature, albedo, downwelling surface short and long-wave fluxes were processed by LandSAF. It was observed that the difference between these parameters varies systematically over the day, together with the dependency on wind speed. Since land surface temperature varies significantly over time and space, a downscaling procedure to a 1000 m spatial resolution was performed. The method has been successfully tested for two areas (Slovenia and southern Germany) during the second half of 2005. The results are slightly better in Slovenia (RMSE = 1.8 °C) than in Germany (2.1 °C). However, both provided a promising result especially considering the high temporal (30 min) and spatial resolution (1000 m) of the results. A “side product” of the thesis is also an LST model of high spatial and temporal resolution that was tested in Slovenia. LST can be measured in a high temporal resolution merely from geostationary satellites which have a poor spatial resolution. Since LST is a spatially inhomogeneous parameter, it was downscaled using contextual analyses and wavelet transform through which we managed to also obtain a high spatial resolution. Already previous studies have shown that the SEVIRI LST is usually greater than the MODIS LST, thus a SEVIRI LST correction was made in order to fit this data to MODIS, because it is operational for longer and presumably more accurate than SEVIRI. The bias was removed after applying the correction; the standard deviation of the corrected LST equals 1.6 °C. Afterwards the downscaling procedure was applied. The first step of downscaling LST from a 5000 m to a 1000 m spatial resolution was performed through contextual analyses – searching for the most appropriate explanatory variables within a moving window (in 5000 m resolution) and applying the multiple regression to the high resolution data (in 1000 m resolution) corresponding to the chosen explanatory variables in order to determine a LST of a higher resolution. The results were tested with the MODIS LST and the average RMSE equalled 2 °C. The second downscaling step was performed by wavelet transform, which was used to fuse the low input frequency of 1000 m LST with the high frequency solar radiance exposure data modelled in a 125 m spatial resolution. The LST computed from the sixth channel of Landsat TM 5 image was used to estimate the quality of the results. It was Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. discovered that the 1000 m spatial resolution provides more reliable data comparing to the 125 m LST, because the downscaling procedure provides good results only in areas with a rough relief and homogeneous land cover. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. POVZETEK Vertikalni energijski tokovi v bližini površja so odvisni od temperature tal. Zato zrak ni ogrevan neposredno od Sonca, ampak posredno od tal. Podobno velja tudi za ohlajanje ob (nočnem) termičnem sevanju površja. Temperatura tal je merjena na nekaterih meteoroloških postajah na različnih globinah, vendar se te meritve ne izvajajo na vseh postajah in tudi ne na sami površini tal. Je prostorsko dokaj nehomogena meteorološka spremenljivka, ki ima značilen povprečni dnevni in letni hod. Najvišje dnevne vrednosti doseže v povprečju zgodaj popoldne in najnižje tik pred sončnim vzhodom. Temperaturo površja (tanka gornja plast rastja), ki se lahko razlikuje od temperature tal (prst ali kamnina), lahko določimo tudi na podlagi daljinskega zaznavanja. Najprimernejši satelitski senzorji za določanje temperatur površja v srednji in visoki prostorski ločljivosti so SEVIRI (MSG), MODIS (Terra in Aqua), AVHRR (NOAA) in ASTER (Terra). Ti senzorji delujejo v termičnem (srednjem in dolgovalovnem) infrardečem spektru elektro-magnetnega valovanja (EMV) in na podlagi tako zajetih podob lahko določimo temperaturo površja (Wan in Dozier, 1996). Prostorska ločljivost omenjenih podob na nehomogenih območjih ne zadošča potrebam nekaterih aplikacij, vendar je bilo do sedaj le malo raziskav usmerjenih v izboljšanje prostorske ločljivosti temperature površja, pridobljene z daljinskim zaznavanjem (Venkateshwarlu in drugi, 2004). Podobno kot temperatura tal doseže temperatura zraka najvišje dnevne vrednosti zgodaj popoldne in najnižje tik pred sončnim vzhodom. Je linearno odvisna od nadmorske višine: v povprečju pade 6,5 °C na vsakih 1000 m v prosti atmosferi (Rakovec in drugi, 1998), nad tlemi na območju Slovenije pa nekako do 4 °C na 1000 m (Mekinda-Masaron, 1995). Lokalni Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. vpliv na temperaturo je najbolj izrazit pri najnižjih temperaturah. Pri veliki prostorski variabilnosti je obstoječ sistem meteoroloških postaj v Sloveniji za določene aplikacije (npr. za ugotavljanje območij zmrzali) preredek, saj je značilna horizontalna dimenzija variabilnosti temperaturnega polja (ponekod manj kot kilometer) v tem primeru precej manjša od povprečne razdalje med meteorološkimi postajami (nekaj deset kilometrov). V primerih, ko je lokalni vpliv nanjo močan (omejeno mešanje zraka ob nizkih temperaturah, lokalno pregrevanje zraka ob močnem sončnem obsevanju), bi morale biti temperature merjene časovno in prostorsko dovolj pogosto, da bi zadostile kriterijem uporabe v aplikacijah srednje in visoke prostorske ločljivosti. Vzpostavitev, vzdrževanje in tudi samo nadzorovanje zgoščenega sistema merilnih postaj zahteva veliko sredstev, zato je bolj smiselno modelirati temperature na podlagi drugih dejavnikov v prostoru, ki vplivajo na njihovo prostorsko porazdelitev. Gre za modeliranje prostorskih spremenljivk, zato je najprimerneje uporabiti metode geografskih informacijskih sistemov (Dolinar, 2004). V Sloveniji so bila npr. izračunana polja za dolgoletne povprečne mesečne in letne temperature (obdobje 1961–1990) s prostorsko ločljivostjo 100 m (Marolt, 2000). Temperaturna polja za posamezne situacije ali povprečne mesečne temperature se računajo s kombiniranim pristopom (geostatistični in fizikalni model), navadno v prostorski ločljivosti približno 1 km (Dolinar, 2004). Do sedaj je bilo v svetu razvitih nekaj diagnostičnih modelov, ki določajo temperature zraka v visoki prostorski ločljivosti na podlagi zgoščene mreže meteoroloških opazovalnic (Brossard in drugi, 2002; Joly in drugi, 2003; Pape in Löffler, 2004). Da določimo temperaturo zraka na podlagi temperature površja, moramo imeti na voljo še dodatne podatke in razumeti povezavo med sevalno temperaturo površja in temperaturo zraka (Gallo in drugi, 1993; Kawashima in drugi, 2000; Garand in drugi, 2004; Meteotest, 2004; Sun in drugi, 2005). 8.1 Hipoteza Cilj naloge je izdelava metodologije, s katero je možno izdelati zvezno polje temperature zraka v visoki prostorski in časovni ločljivosti. To polje je mogoče pridobiti z interpolacijo meritev ali pa s parametrizacijo temperature iz neodvisnih podatkov. Interpolacija je možna, Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. kadar je na voljo dovolj meritev. V Sloveniji so meritve temperature zraka v okviru naloge na voljo le ob 7:00, 14:00 in 21:00, zato iz njih ne moremo dobiti več kot tri zvezna polja na dan. Zato je bila temperatura zraka parametrizirna iz temperature površja in drugih podatkov, pridobljenih z daljinskim zaznavanjem. Temperatura površja se lahko meri z visoko časovno ločljivostjo le iz geostacionarnih satelitov, ki pa imajo slabo prostorsko ločljivost. Ker je temperatura površja prostorsko nehomogena spremenljivka, je v nalogi najprej izboljšana njena ločljivost s kontekstualnimi analizami in valčno transformacijo ob uporabi podatkov boljše prostorske ločljivosti. V skladu s predstavljenimi možnostmi modeliranja so bile v nalogi testirane štiri hipoteze. Prva hipoteza: »Na podlagi podob geostacionarnih satelitov, ki podajajo temperaturo površja v srednji prostorski in visoki časovni ločljivosti (SEVIRI 5000 m in 15 min), je možno modelirati temperaturo površja v veliki prostorski (podobno kot na spletu dosegljivi produkti temperature površja, torej 1000 m ali boljše) in časovni ločljivosti (15 min) z uporabo ustreznih podatkov.« Druga hipoteza: »Temperaturo zraka na višini dveh metrov v srednji prostorski ločljivosti lahko določimo iz temperature zemeljskega površja, če imamo na voljo dodatne podatke (albedo, dolgovalnovno sevanje itd.), ki pojasnjujejo razlike med temperaturo površja in zraka. Postopek omogoča boljšo točnost rezultatov, kot je točnost vremenske napovedi, ki znaša 2 do 3 °C (ARSO, 2007b; DWD, 2007).« Tretja hipoteza: »V nalogi predstavljen model prenosa temperature površja v temperaturo zraka temelji le na konvekciji (vertikalno gibanje zračnih mas), advekcija (horizontalno gibanje zračnih mas) je zanemarjena, ker v okviru naloge ni bilo na voljo ustrezni vetrnih polj. Zato je točnost postopka določitve temperature zraka iz temperature površja bolj zanesljiva v mirnem brezvetrnem vremenu.« Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Četrta hipoteza: »Ko so na voljo meritve temperature zraka, jih je možno interpolirati v zvezno polje visoke prostorske ločljivosti z upoštevanjem statičnih podatkov, ki opisujejo relief in rastje. Postopek omogoča boljšo točnost rezultatov, kot jo omogoča vremenska napoved, ki znaša 2 do 3 °C (ARSO, 2007b; DWD, 2007).« Hipoteze so bile preizkušene po različnih poglavjih. Po uvodu je v drugem poglavju naloge preizkušena četrta hipoteza. V tem poglavju je tudi nakazano, da omogoča velika prostorska ločljivost načeloma boljše rezultate kot majhna ločljivost, zato je v tretjem poglavju opisan postopek izboljšave ločljivosti s kontekstualnimi analizami s 5000 m na 1000 m in v četrtem še dodatna izboljšava z valčno transformacijo s 1000 m na 125 m (prva hipoteza). V petem poglavju je nato določena temperatura zraka na podlagi temperature površja najprej prostorski ločljivosti 5000 m in nato še v izboljšani ločljivosti 1000 m. S tem je preizkušena druga hipoteza. Na koncu je v zaključku analiziran vpliv vetra na odstopanja, s čimer je preskušena tretja hipoteza. 8.2 Interpolacija temperature zraka na nehomogenih območjih Mnogo študij potrjuje dejstvo, da imata relief in rastje velik vpliv na prostorsko porazdelitev temperature površja (Anquetin et al., 1999; Benichou and Le Breton, 1987; Bolstad et al., 1998; Bootsma, 1976; Joly et al., 2003; Geiger, 1965; Yoshino, 1975). Največji vpliv ima prav gotovo nadmorska višina (temperatura zraka pada z nadmorsko višino zaradi manjšega zračnega pritiska). Te vplive lahko upoštevamo pri interpolaciji meritev v zvezno polje. Uporabimo lahko več različnih metod, ki so nam na voljo v raznih GIS programskih paketih. Interpolacijske metode se v grobem delijo na deterministične in geostatistične (Kanevski and Maignan, 2004). Obe vrsti temeljita na visoki stopnji prostorske avtokorelacije – predvideva se, da so si območja, ki so si bližje, bolj podobna, kot tista, ki so med sabo bolj oddaljena. Deterministične metode so hitre, geostatistične počasnejše, a minimizirajo šum v podatkih. Kadar je prostorska avtokorelacija nizka, je treba k izdelavi zveznega temperaturnega polja Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. pristopiti drugače: prostorsko porazdelitev temperature zraka je mogoče določiti tudi iz nekaterih naravnih parametrov, ki jih lahko razumemo kot atribute merilnih postaj. Če je povezava med njimi in temperaturo zraka statistično pomembna, potem vplivajo atributi na prostorsko porazdelitev temperature zraka. Zato določimo zvezo med vsemi statistično pomembnimi atributi in temperaturo zraka na merilnih postajah in jo nato uporabimo s podatki, zapisanimi v pravilni mreži. Na ta način (multipla regresija) določimo prostorsko porazdelitev temperature za vse območje. Tak pristop ima kar izbrane prednosti pred klasičnimi oblikami interpolacije, saj lahko na podlagi znanja o vplivih na temperaturo zraka predvidimo atribute, ki jih nato uporabimo le, kadar so statistično povezani s temperaturo. Sicer pa je lahko ta pristop zelo počasen, če imamo na voljo mnogo meritev in mnogo pojasnjevalnih atributov. Uporabljena metoda je sestavljena iz sedmih korakov (slika 2): • izračun prostorske avtokorelacije meritev temperature zraka (če je različna od nič, uporabimo deterministične ali geostatistične metode, sicer nadaljujemo), • izbira pojasnjevalnih atributov, • priprava podatkov, • multipla regresija, • ocena rezultatov z navzkrižnim preverjanjem, • dodatna interpolacija odstopanj in • kartiranje rezultatov. Metoda je bila preskušena v Sloveniji, kjer so bile v okviru naloge na voljo meritve temperature zraka na višini dveh metrov na dvajsetih postaj za leto 2005 ob 7:00, 14:00 in 21:00 po srednjeevropskem času. Njihova prostorska porazdelitev je bila pojasnjena s podatki, izpeljanimi iz digitalnega modela višin DMR 100 (Podobnikar & Mlinar, 2006), podatki o vrsti rabe tal (Kokalj in Oštir, 2006) in MODIS-ovimi satelitskimi podatki: Podatki izpeljani iz DMV Podatki izpeljani iz vrste rabe tal MODIS-ovi satelitski podatki višine umetno izdelan vegetacijski indeks NDVI nakloni reliefa razdalje do morja EVI usmerjenosti reliefa razdalje do gozda temperatura površja kvaziglobalna osončenost razgibanost reliefa TduP (višinska razlika od okolice) Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Ker imajo nekateri od teh atributov lokalni in nekateri regionalni vpliv, so bili mnogi od njih določeni znotraj okolic različnih velikosti (preglednica 1). V prvem poskusu (Test 1) MODIS-ovi satelitski podatki in razdalja do morja še niso bili na voljo. Od preostalih podatkov so bili vsi (razen razdalje do gozda in TduP) določeni znotraj okolic sedem različnih velikosti (preglednica 1), tako da je bilo na koncu na voljo 44 slojev primernih za pojasnitev prostorske porazdelitve temperature. Pokazalo se je, da se vpliv določenih atributov zelo spreminja v odvisnosti od velikosti območja, v katerem je bil atribut določen (slika 7). Rezultati Testa 1 niso bili zadovoljivi zaradi nekaj zelo velikih odstopanj (sicer je standardni odklon za vse leto znašal 1,8 °C), zato se je bilo treba vrniti nekaj korakov nazaj – k izbiri pojasnjevalnih atributov. Ker so bili največja odstopanja na postajah blizu morju, je bil kot dodaten atribut izbrana oddaljenost do morja. Skupna točnost se je popravila za desetinko stopinje, velika odstopanja na postaji letališče Portorož pa so ostala. Ugotovljeno je bilo, da gre za napačen položaj postaje – glede na podatke vrste rabe tal se naj bi postaja nahajala v morju, zato je bila »premaknjena na suho«. Poleg tega pri izračunu oddaljenosti do gozda najprej ni bilo upoštevano grmičevje, zato je bil izveden poskus tudi z upoštevanjem grmičevja kot gozda in še brez upoštevanja oddaljenosti do gozda. Ugotovljeno je bilo, da oddaljenost do gozda nima velikega vpliva na porazdelitev temperature, zato je bil atribut izpuščen iz analize in s tem se je skupna točnost povečala za 0,05 °C. Končna točnost 1,6 °C je bila dosežena z dodatnim upoštevanjem MODIS-ovih satelitskih podatkov (NDVI, EVI in temperature površja). Slika 13 prikazuje primer končnih rezultatov. Rezultati so primerljivi s predhodnimi študijami (npr. Anderson, 2002), ne pa tudi npr. s študijo Jolya in sodelavcev (2003), ki so dosegli natančnost 0,3 °C, vendar le na območju velikosti 8 km2 ob 50 merilnih postajah. Rezultati sicer niso bili preverjeni s popolnoma neodvisnimi podatki, ker je bilo merilnih postaj pač premalo. Za oceno točnosti je bila zato uporabljena metoda navzkrižnega preverjanja, kjer je v okviru vsake interpolacije ena od meritev izpuščena iz izračuna. Postopek izračuna odstopanja se ponovi tolikokrat, kolikšno je število postaj. Na podlagi vseh odstopanj se nato oceni točnost rezultatov. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 8.3 Izboljšava prostorske ločljivosti z uporabo kontestualnih analiz V drugem poglavju je bilo ugotovljeno, da je temperatura površja eden od tistih parametrov, ki pomembno vplivajo na temperaturo zraka. Poleg tega je tudi jasno, da je prostorska ločljivost temperature površja, ki jo posname SEVIRI (na območju Slovenije znaša približno 5000 m), preslaba za kakovostno določitev temperature zraka. Zato je bila v tretjem poglavju ločljivost temperature površja izboljšana s kontestualnimi analizami ob uporabi podatkov boljše prostorske ločljivosti. Na področju izboljšave ločljivosti površja še ni bilo predstavljenih veliko študij, razen kadar gre za izboljšavo multispektranih s pankromatskimi posnetki, na področju izboljšave prostorske ločljivosti temperature površja pa (razen Venkateshwarlu in drugi, 2004) literature ni. Primerjava SEVIRI-jeve temperature površja z ostalimi primerljivi produkti je pokazala, da je SEVIRI-jeva temperatura v povprečju nekaj stopinj višja kot npr. MODIS-ova temperatura površja (Madeira in drugi, 2005; Noyes in drugi, 2006). Razlike bi bile lahko posledica v različni karti emisivnosti, ki je uporabljena za izračun ene ali druge temperature, glavni razlog pa je najbolj verjetno kot gledanja – SEVIRI je na krovu geostacionarnega MSG, ki gleda na Slovenijo vedno z juga, zato vidi več proti jugu obrnjenih pobočij in manj senc kot MODIS, ki preleti območje skoraj v nadirju (slika 15). Zato je bila najprej izvedena korekcija SEVIRIjeve temperature v izvirni ločljivosti. Poznavanje omenjenih študij in vizualna primerjava z ostalimi podatki sta pokazala, da so razlike odvisne od: • deklinacije in zenitnega kota sonca, • albeda, • deleža pokritosti z rastjem (FVC), • nadmorske višine in • popravkov glede na obliko reliefa ter rabe tal površja. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Z uporabo multiple regresije je bil na podlagi 14 posnetkov MODIS-a (preglednica 6) izdelan popravek za SEVIRI-jevo temperaturo površja (ocenjena natančnost znaša 1.6 °C). Popravljena temperatura površja je bila potem končno pripravljena na izboljšanje ločljivosti. Kontekstualne analize so se izkazale kot primerno orodje za izboljšavo ločljivosti. V preteklosti se je namreč že izkazalo, da je temperatura površja linearno povezana z vegetacijskim indeksom NDVI (slika 20). Povezava ni konstanta, saj se spreminja v odvisnosti od časa in kraja. Zato je najbolj smiselno uporabiti kontekstualno analizo, pri kateri določimo parametre regresijske premice med temperaturo površja in NDVI znotraj premikajočega se okna, poljubne velikosti. Regresijski parametri so torej določeni v ločljivosti 5000 m in če jih prevzorčimo na ločljivost 1000, lahko z uporabo NDVI v prostorski ločljivosti 1000 m izračunamo temperaturo površja v isti ločljivosti (slika 21). Ker je temperatura površja odvisna tudi od drugih parametrov, kot so nadmorska višina, albedo, vpadni kot sonca in EVI, lahko uporabimo multiplo regresijo, kjer iz statistično pomembnih parametrov za vsako rastrsko celico velikosti 5 krat 5 km določimo svojo enačbo. Največja prednost te metode je preprosta uporaba, slabost pa počasno delovanje. Izkazalo se je, da znaša korelacijski indeks med izboljšano temperaturo površja in MODIS-ovo temperaturo običajno 0,8. Pričakovana natančnost je enaka SEVIRI-jevi natančnosti (povprečni RMSE znaša 2 °C). Slika 30 prikazuje primer, kako se rezultati ujemajo z referenčnimi MODISovimi podatki – vizualna primerjava pokaže, da sta si obe temperaturi v relativnem smislu dokaj podobni (topla območja so topla na obeh slikah), statistična primerjava pa vseeno pokaže razlike v absolutnem smislu, saj so na novo izdelane temperature bolj generalizirane. Razlike so tudi posledica nepopolnega ujemanja temperatur v prostorski ločljivosti 5000 m. 8.4 Izboljšava prostorske ločljivosti z uporabo valčne transformacije Ena ob možnih metod za izboljšavo prostorske ločljivosti je tudi združevanje podob z valčno transformacijo, ki je bila že večkrat uporabljena za izboljšavo multispektralnih s pankromatskimi satelitskimi podobami (Ranchin in drugi, 2003; Acerbi-Junior in drugi, 2006). Vendar temeljijo vse pretekle študije na združevanju podobnih si podatkov, v nalogi pa Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. je preizkušeno, kako lahko z valčno transformacijo združimo temperaturo površja v prostorski ločljivosti 1000 m s podatki o simulirani osončenosti v prostorski ločljivosti 125 m. Izboljšava temperature površja z osončenostjo se zdi logična izbira v primeru razgibanega reliefa, na ravninah pa na prostorsko porazdelitev temperature vplivajo drugi dejavniki, kot je npr. vlažnost tal. Zato je bilo za testno območje izbrano območje v Alpah na meji med Slovenijo in Avstrijo, velikosti 32 krat 32 km. Na osončenost vplivajo različni dejavniki: položaj sonca na nebu, vremenske razmere in vpliv oblike reliefa. Ker so na testnem območju vremenske razmere dokaj podobne, položaj sonca pa praktično enak, je bilo predvidevano, da na osončenost vpliva le oblika reliefa, zato so bili na podlagi digitalnega modela (Podobnikar in Mlinar, 2006) višin pripravljeni naklon ter usmerjenost reliefa in delež vidnega neba. Naklon in usmerjenost reliefa sta pomembna za določitev vpadnega kota sonca, od katerega je odvisna direktna osončenost na tleh (slika 34). Če je vpadni kot sonca večji od 90°, je območje v lastni senci. Območje je lahko tudi v vrženi senci, kadar je med območjem in soncem ovira. Ker območja v senci prejmejo le difuzno osončenost, je zelo pomembno določiti, katera območja so obsijana s soncem in katera ne. Oblika reliefa vpliva tudi na difuzno osončenost. Območja na ravnini vidijo namreč mnogo več svetlega neba kot območja v zaprtih dolinah. Delež difuzne osončenosti, ki jo prejmejo območja v dolinah, je zato premo sorazmeren deležu vidnega neba. Del difuzne osončenosti pride tudi od bolj ali manj svetlih tal, zato je treba pri izračunu upoštevati še albedo površja. Na podlagi klimatskih podatkov (Kastelec in drugi, 2007) je bilo ocenjeno razmerje med difuzno in direktno osončenostjo, kar je omogočilo, da se oceni vpliv oblike reliefa kot vrednost, ki je premo sorazmerna dejanski osončenosti. Ta vrednost je bila nato z linearnim raztegom histograma prilagojena histogramu temperature površja, kar je omogočilo združevanje iz osončenosti določenega približka temperature v prostorski ločljivosti 125 m s temperaturo površja v prostorski ločljivosti 1000 m (slika 39). Z diskretno valčno transformacijo so bili izračunani valčni koeficienti za oba podatka. Izkazalo se je, da so rezultati najboljši, če pri valčni transformaciji uporabimo Daubechie-jevo družino vlačkov pete stopnje. Po transformaciji (do četrte stopnje) so bili koeficienti obeh podatkov na različnih stopnjah povprečeni: Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. • na četrti stopnji so bili kot izhodni koeficienti uporabljeni le koeficienti iz kilometrske temperature površja, • na naslednjih stopnjah pa so bile uteži: 0,7 in 0,3 (3.), 0,5 in 0,5 (2.), 0,3 in 0,7 (1. stopnja); prva utež v paru velja za temperaturo v ločljivosti 1000 m in druga za temperaturo ločljivosti 125 m. Te uteži so dale najboljši rezultat – najboljše ujemanje z referenčno temperaturo dne 30. 7. 2005 ob 10.30. Takrat je namreč testno območje posnel senzor Landsat TM 5, ki je zajel podatke s senzorjem TM, in na podlagi podatkov šestega kanala omenjenega senzorja je bila izračuna temperatura površja v prostorski ločljivosti 120 m (Oštir, 2006). Statistična primerjava rezultatov z Landsat-ovo temperaturo je pokazala dokaj velika odstopanja, ki znašajo ±7 °C, korelacijski indeks 0,71, RMSE pa 1,5 °C. Vizualna analiza (slika 40) je sicer pokazala, da so po pričakovanjih soncu bolj izpostavljena območja toplejša od tistih v senci, po drugi strani se je tudi izkazalo, da je poprečna temperatura po izboljšanju ločljivosti na 125 m previsoka glede na referenčne podatke. Razlike so tudi posledica neupoštevanja različne rabe tal, ki ima na ravnih območjih zelo velik vpliv na temperaturo površja. Največji praktični problem v postopku je prilagoditev histograma osončenosti na histogram temperature površja v kilometrski ločljivosti. Ta dva podatka sta si preveč različna, zato bi bilo treba modelirati temperaturo površja (ob upoštevanju osončenosti ne bi smeli zanemariti fizičnih lastnosti površja). Glede na vse težave je bilo na koncu odločeno, da preizkušene metode nima smisla uporabljati, ker je vpliv drugih dejavnikov prevelik, da bi lahko samo s simulirano osončenostjo izboljšali ločljivost temperature površja. 8.5 Parametrizacija temperature zraka iz temperature površja Ker se zrak segreva posredno od tal, je temperatura površja, zaznana s sateliti, primeren podatek za določitev temperature zraka. Analiza vremenske napovedi temperature zraka v Sloveniji za 24 ur naprej je pokazala, da znaša točnost napovedi približno 2–3 °C (najvišje in najnižje temperature so načeloma bolje določene), težave pa predstavljajo predvsem nekatere Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. posebne vremenske situacije, kot je npr. megla (takrat je lahko točnost tudi trikrat slabša). Podobne rezultate dosegajo tudi v Nemčiji (slika 42). Določitev temperature zraka na podlagi temperature površja je ena še nerešenih nalog sodobne meteorologije. Tudi na tem področju ni veliko objav kljub temu, da so se prve študije na to temo pojavile že pred skoraj tridesetimi leti. Objavljene metode lahko razdelimo v štiri skupine. V prvo skupino sodijo tista dela, kjer je temperatura zraka le linearno povezana s temperaturo površja (Chen in drugi, 1983; Davids in Tarpley, 1983; Greena in Hay, 2002). Te študije dajejo sicer dobre rezultate, vendar le v nekaterih posebnih okoliščinah). Tudi druga skupina študij ima statistično osnovo, vendar ne gre več le za osnovno regresijo, saj upoštevajo tudi npr. zenitni kot sonca (Cresswell, 1999), Jang (2004) je delal z nevronskimi mrežami itd. Tudi TVX metoda je v osnovi statistična metoda, vendar je bila objavljena in preizkušena večkrat (Saravanapavan and Dye 1995, Prihodko and Goward 1997, Czajkowski et al. 2004; Zidarič, 2007). Temelji na dejstvu, da več rastja načeloma pomeni nižjo temperaturo površja. Poleg tega lahko še predvidevamo, da znaša temperatura znotraj neskončno debele krošnje natanko toliko kot temperatura vrha krošnje – na kratko, temperatura zraka je kar enaka s satelitom zaznani temperaturi površja v primeru veliko vegetacije. Če ocenimo parametre regresijske premice na osnovi kontekstualne analize (slika 43), lahko ob izbrani vrednosti NDVI-ja »neskončno debele krošnje« izračunamo temperaturo zraka, ki običajno dosega točnost 3–4 °C. Četrta skupina študij ima fizikalno ozadje (Meteotest, 2004; Pape in Löffler, 2004; Sun in ostali, 2005). Ti pristopi upoštevajo energetsko ravnotežje na zemeljskem površju z modeliranjem energijskih tokov (slika 44). Načeloma so bolj točni od ostalih modelov (največja odstopanja so manjša od 5 °C), vendar pa za svoje delovanje potrebujejo veliko podatkov. V nalogi sta preizkušena dva od predstavljenih pristopov: • TVX zaradi svoje preprostosti in • Meteotest-ova metoda, ki je edina operativna metoda (znotraj programskega paketa Meteonorm). Točnost TVX pristopa dosega na območju Slovenije 2,5 °C pri uporabi MODIS-ove temperature površja in desetinko stopinje slabše rezultate pri SEVIRI-jevi temperaturi površja Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. v izvirni prostorski ločljivosti. Uporaba Meteotest-ove parametrizacije je sprva prinesla slabe rezultate, kajti koeficienti parametrizacije so bili kalibrirani le na petih švicarskih postajah. Po kalibraciji parametrov s slovenskimi meritvami se je RMSE znižal na 2,3 °C. Rezultati so bili mnogo slabši v Nemčiji, kjer je RMSE znašal na 3,8 °C. Ugotovljeno je bilo, da v nekaterih primerih prostorska ločljivost podatkov omejuje točnost rezultatov (odstopanja na postajah Golica v Sloveniji in Zugspitze v Nemčiji). Zato je bila ločljivost izboljšana na osnovi kontekstualne analize (najprej samo z NDVI-jem). Poleg tega so bila odstopanja odvisna od časa – največja okoli poldneva in najnižja pred oziroma po sončnem zahodu. V ta namen je bila izdelana nova parametrizacija, ki je temeljila na razlikah med temperaturo površja in temperaturo zraka. Te so bile pojasnjene z multiplo regresijo naslednjih spremenljivk: • kratkovalovnega sevanja (popravljenega za vpliv albeda), • dolgovalovnega sevanja, • kosinusa zenitne razdalje množenega z naravnim logaritmom NDVI in • vetra. Točnost v Sloveniji (RMSE = 1,8 °C; slika 65) in Nemčiji (RMSE = 2,1 °C) se je z novo parametrizacijo močno izboljšala, vendar se koeficienti za obe območji razlikujejo. Preizkušena je bila tudi parameterizacija, osnovana na izboljšani temperaturi tal v ločljivosti 1000 m (opisano v 3. poglavju) brez vetra, a se točnost ni izboljšala, ampak poslabšala. Možen razlog za to je popravek za SEVIRI-jevo temperaturo površja, ki je bil izdelan v večini na osnovi dnevnih posnetkov MODIS-a na satelitu Aqua, ki naše kraje prečka okoli 12:00 GMT. Zato je bila parameterizacija preizkušena še enkrat ob 13:00 GMT in rezultati znašajo: minimum -3,7 °C, maksimum 5,4 °C, RMSE = 1,8 °C, r = 0,97. Izpopolnjen model temperature tal je torej smiselno uporabiti podnevi, medtem ko v splošnem lahko dosegamo dobre rezultate tudi z izboljšavo prostorske ločljivosti samo na osnovi NDVI. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 8.6 Zaključek Za konec sta bili med sabo primerjani še obe opisani metodi določitve temperature zraka. Najpomembnejša je ocena točnosti, kjer lahko ugotovimo, da zagotavlja interpolacija boljše rezultate na večini postaj po Sloveniji (slika 69). Poleg tega je bilo preverjeno, kolikšen vpliv ima hitrost vetra na določitev temperature zraka. Ugotovljeno je bilo, da je točnost rezultatov pri interpolaciji neodvisna od hitrosti vetra, pri prametrizaciji iz temperature površja pa je točnost metode manjša pri močnejšem vetru. Modeliranje temperatur površja in zraka na večjih območjih je bilo do sedaj izvedeno le v slabi prostorski ali časovni ločljivosti. Rezultat naloge je metodologija, s katero lahko v geografskih informacijskih sistemih modeliramo temperature površja v srednjem in velikem merilu v veliki časovni ločljivosti na podlagi podatkov daljinskega zaznavanja, podatkov o reliefu in rastju. Model temperature površja v prostorski in časovni ločljivosti lahko v nekaterih vremenskih situacijah (kadar je veter šibak in površje ni pokrito s snegom) uspešno uporabimo podnevi tudi za določitev temperature zraka. Zvezno polje temperature zraka je bilo sicer določeno tudi z interpolacijo, s katero so bili rezultati boljši, vendar je interpolacija mogoča le ob zadostnem številu meritev. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Literature and sources Acerbi-Junior, F.W., Clevers, J.G.P.W., Schaepman M.E. 2006. The assessment of multi-sensor image fusion using wavelet transforms for mapping the Brazilian savanna. International journal of applied earth observation and geoinformation 8, 4: 278–288. Agnew, M.D., Palutikof, J. 1996. GIS-based downscaling of climate data for the Mediterranean basin using terrain Information. In Bindi, M., Gozzini, B. (ed.) Seminar on data spatial distribution in meteorology and climatology. Luxembourg, Office for official publications of the European community: 201– 206. Allen, R.G., Trezza, R., Tasumi, M. 2006. Analytical integrated functions for daily solar radiation on slopes. Agricultural and forest meteorology 139, 1–2: 55–73. Anderson, S. 2002. An evaluation of spatial interpolation methods on air temperature in Phoenix, AZ. http://www.cobblestoneconcepts.com/ucgis2summer/anderson/anderson.htm (12.12.2006). Anquetin S, Guilbaud C, Chollet JP, 1999. Thermal valley inversion on the dispersion of a passive polluant in a complex mountainous area. Revue atmospheric environment 33: 3953–3959. Ansari, A.Q., Loomis, W.E. 1959. Leaf temperatures. American journal of botany 46: 713–717. ARSO 2007a. ALADIN/SI – napoved vetra z ločljivostjo 2,5 km. http://www.arso.gov.si/vreme/napovedi%20in%20podatki/dada/. (12.3.2007). ARSO 2007b. Članek o verifikaciji regionalnih vremenskih napovedi na ARSO. http://www.arso.gov.si/podro~cja/vreme_in_podnebje/poro~cila_in_projekti/Cl anek_verifikacija.pdf (12.3.2007). Baldocchi, D. 2006. Lecture 8, Solar radiation, Part 2, Earth-Sun geometry. http://nature.berkeley.edu/biometlab/espm129/pdf/lecture%209%20espm%201 29.pdf (12.3.2007). Benichou, P., Le Breton, O., 1987. Prise en compte de la topographie pour la cartographie des champs pluviométriques statistiques. La météorologie 7, 19: 22–34. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Bolstad, P., Swift L., Collins F., Regniere J., 1998. Measured and predicted air temperatures at basin to regional scales in the southern Appalachian mountains. Agricultural and forest meteorology 91: 161–176. Bootsma, A., 1976. Estimating minimum temperature and climatological freeze risk in hilly terrain. Agricultural meteorology 16: 425–443. Brossard, T., Elvebakk, A., Joly, D., Nilsen, L. 2002. Modelling index of thermophily by means of a multi-source database on Broggerhalvoya Peninsula (Svalbard). International journal of remote sensing 23: 4683–4698. Burrough, P., McDonnell, R. A. 1998. Principles of geographical information systems. Oxford, Oxford university press: 327 p. Chapman, L., Thornes, J.E. 2003. The use of geographical information system in climatology and meteorology. Progress in physical geography 27: 313–330. Chen, E., Allen, L.H.Jr., Bartholic, J.F., Gerber, J.F. 1983. Comparison of winter- nocturnal geostationary satellite infrared-surface temperature with shelter- height temperature in Florida. Remote sensing of environment 13, 4: 313–332. Chen, S., Zhu, H.Y. 2007. Wavelet transform for processing power quality disturbances. EURASIP journal on advances in signal processing: 1–20. Chung, U., Yun, J.I. 2004. Solar irradiance-corrected spatial interpolation of hourly temperature in complex terrain. Agricultural and forest meteorology 126: 129– 139. Collins, F.C., Bolstad, P.V. 1996. A comparison of spatial interpolation techniques in temperature estimation. In Proceedings of the third international conference / workshop on integrating GIS and environmental modeling CD-ROM (Santa Fe, 1996). http://www.ncgia.ucsb.edu/conf/SANTA_FE_CDROM/ sf_papers/collins_fred/collins.html (12.12.2006). COST 2000. Draft memorandum of understanding; the use of geographic information system in climatology and meteorology. http://www.knmi.nl/samenw/cost719/documents/mou.pdf (12.12.2006). Cresswell, M.P. 1999. Estimating surface air temperatures from METEOSAT land surface temperatures using an empirical solar zenith angle. International journal of remote sensing 20, 6: 1125–1132. Czajkowski, K.P., Goward, S.N., Mulhern, T., Goetz, S.J., Walz, A., Shirey, D., Stadler, S., Prince, S.D., Dubayah, R.O. 2000. In Dale, A., Luvall, J.C. (ed.) Estimating environmental variables using thermal remote sensing. Thermal remote sensing in land surface processes. New York, CRC press: 11–32. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. DAAC, 2006. EOS data gateway. http://edcimswww.cr.usgs.gov/pub/imswelcome/ (12.12.2006). Dash, P. 2004. Land surface temperature and emissivity retrieval from satellite measurements. PhD thesis. Univeristy of Karlsruhe: 77 f. Davis, F.A., Tarpley, J.D. 1983. Estimation of shelter temperatures from operational satellite sounder data. Journal of applied meteorology 22, 3: 369–376. Dobesch, H., Tveito, O.E., Bessemoulin, P. 2001. Geographical information system in climatological application, KLIMA report 13/01. Oslo, Norwegian meteorological institute: 67 f. Dolinar, M. 2004. GIS kot orodje pri izdelavi klimatskih kart. In Podobnikar, T., Perko, D., Hladnik, D. Krevs, M., Čeh, M., Stančič, Z. (ed.) Geografski informacijski sistemi v Sloveniji 2003–2004. Ljubljana, Založba ZRC: 195– 202. DWD 2007. Verifikation numerischer Wettervorhersagen. http://www.dwd.de/de/FundE/Analyse/Verifikation/Verifikation.htm (12.3.2007) EARS 2003. Air temperature mapping. http://www.ears.nl/ (12.3.2007). EEA 2005. Corine land cover 2000. http://dataservice.eea.europa.eu/dataservice/metadetails.asp?id=822 (25.3.2007). Efroymson, M.A. 1960. Multiple regression analysis. In Ralston, A., Wilf, H.S. (ed.): Mathematical methods for digital computers. New York, Wiley: 191–2003. Gallo, K.P., McNab, A.L., Karl, T.R., Brown, J.F., Hood, J.J., Tarpley, J.D. 1993. The use of NOAA AVHRR data for assessment of the urban heat island effect. Journal of applied meteorology 32, 5: 899–908. Garand, L., Buehner, M., Wagneur, N. 2004. Background error correlation between surface skin and air temperatures: estimation and impact on the assimilation of infrared window radiances. Journal of applied meteorology 43, 12: 1853–1863. Gardner, T.W., Sasowski, K.C., Day, R.L. 1990. Automated extraction of geomorphometric properties from digital elevation data. Zeitschrift für Geomorphologie 80: 57–68. Geiger, R. 1965. The climate near the ground. Cambridge, Harvard University press: 626 p. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Greena, R.M., Hay, S.I. 2002. The potential of pathfinder AVHRR data for providing surrogate climatic variables across Africa and Europe for epidemiological applications. Remote sensing of environment 79, 2–3: 166–175. Hastings, D.A., Dunbar, P.K. 1998. Development and assessment of the global land one-km base elevation digital elevation model (GLOBE). ISPRS archive 32, 4: 218–221. Heuvelink, G.B.M. 1998. Error propagation in environmental modelling with GIS. New York, Taylor & Frances: 127 p. Hudson, G., Wackernagel, H. 1994. Mapping temperature using kriging with external drift: theory and an example from Scotland. International journal of climatology, 14: 77–91. Hypergeo 2004. Autocorrelation. http://www.hypergeo.eu/article.php3?id_article=223 (12.3.2007). Jang, J.-D. 2004. Evaluation of thermal-water stress of forest in southern Québec from satellite images. PhD thesis. Laval, Faculté de forestrie et de géomatique université Laval Québec: 237 f. Jin, M., Dickinson, R.E. 1999.Interpolation of surface radiative temperature measured from polar orbiting satellites to a diurnal cycle. Journal of geophysical research, 104, 2: 2105–2116. Joly, D., Nilsen, L., Fury, R., Elvebakk, A., Brossard, T. 2003. Temperature interpolation at a large scale: test on a small area in Svalbard. International journal of climatology, 23, 13: 1637–1654. Kanevski M., Maignan M. 2004. Analysis and modelling of spatial environmental data. EPFL press, Lausanne: 300 p. Kastelec, D., Rakovec, J., Zakšek, K. 2007. Sončna energija v Sloveniji. Ljubljana, Založba ZRC: 136 p. Kawashima, S., Ishida, T. Minomura, M. Miwa, T. 2000. Relations between surface temperature and air temperature on a local scale during winter nights. Journal of applied meteorology 39, 9: 1570–1579. Kokalj Ž, Oštir, K. 2006. Ugotavljanje pokrovnosti Slovenije iz satelitskih posnetkov Landsat. Geografski vestnik 78, 2: 85–95. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Kostopoulou, E., Giannakopoulos, C., Flocas, H., Holt, T., Le Sager, P., Psiloglou, B., Hatzaki, M. 2006. Comparison of ERA-40 reanalysis downscaled temperature and precipitation with observational data over Greece. http://www.cosis.net/abstracts/EGU06/05858/EGU06-J-05858.pdf (12.3.2007). Kriebel, K.T., Gesell, G., Kästner, M., Mannstein, H. 2003. The cloud analysis tool APOLLO: Improvements and validation. International journal of remote sensing, 24: 2389–2408. LandSAF 2007. Land surface temperature. http://landsaf.meteo.pt/ (12.3.2007). Laslett, G.M. 1994. Kriging and splines: an empirical comparison of their predictive performance in some applications. Journal of the American statistical association 89: 391–409. Liang, S.S. 2004. Quantitative remote sensing of land surfaces. New York, Willey: 560 p. Lowry, R. 2007. The Significance of a correlation coefficient. http://faculty.vassar.edu/lowry/ch4apx.html (12.3.2007). Madeira, C., Dash, P., Olesen, F., Trigo, I. 2005. Intercomparison of METEOSAT-8 derived LST with MODIS and AATSR similar products. http://www.eumetsat.int/groups/cps/documents/document/pdf_conf_p46_s4_05 _trigo_v.pdf (12.3.2007). Marolt, D. 2000. Prostorske interpolacije minimalne in maksimalne temperature zraka. Diploma thesis. Ljubljana, Faculty of mathematics and physics, Department of physics, Meteorological curse: 90 f. Matheron. G. 1963. Principles of geostatistics. Economic geology 58: 1246–1268. Mekinda-Masaron, T. 1995. Klimatografija Slovenije 1961–1990: temperatura zraka. Ljubljana, Hidrometerorološki zavod Republike Slovenije: 356 p. Méndez, A.A.J. 2004. Estimate ambient air temperature at regional level using remote sensing techniques. Master thesis. Enschede, ITC Enschede: 86 f. Meteotest 2004. Meteonorm handbook, Part III: Theory Part 2. http://www.meteotest.ch/pdf/am/theory_2.pdf (12.3.2007). Noyes, E.J., Corlett, G.K., Kong, X., Remedios, J.J., Llewellyn-Jones D.T. 2006. The AATSR land surface temperature (LST) PRODUCT: comparison with LST from SEVIRI. http://www.leos.le.ac.uk/group/ejn2/FILES/lsa_saf_workshop_2006_paper_aat sr_lst.pdf (12.3.2007). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Oštir, K. 2006. Daljinsko zaznavanje. Ljubljana, Založba ZRC: 250 p. Pajares, G., Manuel de la Cruz, J. 2004. A wavelet-based image fusion tutorial. Pattern recognition 37: 1855–1872. Pape, R., Löffler, J. 2004. Modelling spatio-temporal near-surface temperature variation in high mountain landscapes. Ecological modelling 178, 3–4: 483– 501. Philip. G.M., Watson. D.F. 1986. Matheronian Geostatistics – Quo Vadis? Mathematical geology 18–19: 93–117. Podobnikar, T., Mlinar, J. 2006: Integriranje podatkov reliefa Slovenije. In Perko, D., Nared, J., Čeh, M., Hladnik, D., Krevs, M., Podobnikar, T., Šumrada, R. (ed.) Geografski informacijski sistemi v Sloveniji 2005–2006. Ljubljana, Založba ZRC: 33–41. Polikar, R. 2001. The wavelet tutorial. http://users.rowan.edu/~polikar/WAVELETS/WTtutorial.html (12.3.2007). Potočnik, H. 1928. Das Problem der Befahrung des Weltraums. Berlin, Richard Carl Schmidt: 186 p. Prihodko, L., Goward, S.N. 1997. Estimation of air temperature from remotely sensed surface observations. Remote sensing of environment 60, 3: 335–346. Rakovec, J., Vrhovec, T. 2000. Osnove meteorologije za naravoslovce in tehnike. Ljubljana, Društvo matematikov, fizikov in astronomov Slovenije: 329 p. Ranchin, T., Aiazzi, B., Alparone, L., Baronti, S., Wald, L. 2003. Image fusion—the ARSIS concept and some successful implementation schemes. ISPRS journal of photogrammetry & remote sensing 58: 4–18. Saravanapavan, T., Dye, D.G. 1995. Satellite estimation of environmental variables by the contextual analysis method: validation in a seasonal tropical environment. http://GISdevelopment.net (12.3.2007). Shepard, D. 1968. A two dimentional interpolation function for irregularly-spaced data. In Proceedings of the 23rd ACM national conference. New York, ACM press: 517–523. Short, N.M., Stuart, L.M.Jr. 1982. The heat capacity mapping mission (HCMM) anthology. Washington, NASA: 264 p. StatSoft 2003. Multiple regression. http://www.statsoft.com/textbook/stmulreg.html#index (12.3.2007). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Stisen, S., Sandholt, I. Fensholt, R. 2004. Meteosat second generation data for assessment of surface moisture status. http://earth.esa.int/workshops/msg_rao_2004/papers/27_stisen.pdf (12.3.2007). Sun, Y.-J., Wang, J.-F., Zhang, R.-H., Gillies, R.R., Xue, Y., Bo, Y.-C. 2005. Air temperature retrieval from remote sensing data based on thermodynamics. Theoretical and applied climatology 80, 1: 37–48. Šumrada, R. 2006. Modelling real estate transactions with UML. Geodetski vestnik 50, 4: 597–608. Švab, A., Oštir, K. 2006. High-resolution image fusion: methods to preserve spectral and spatial resolution. Photogrammetry and remote sensing 72, 5, 565–572. Tveito, O.E., Schoner, W. 2002. Applications of spatial interpolation of climatological an meteorological elements by use of geographical information systems (GIS), KLIMA report 28/02. Oslo, Norwegian meteorological institute: 66 f. Venkateshwarlu, C., Gopal, R.K., Prakash, A. 2004. Neural networks in land surface temperature mapping in urban areas from thermal infrared data. IEEE international 3, 20–24: 1589–1590b. Wan, Z., Dozier, J. 1996. A generalized split-window algorithm for retrieving land- surface temperature from space. IEEE transactions on geoscience and remote sensing 34, 4: 892–905. Yoshino M., 1975. Climate in small area. Tokyo: Tokyo university press: 549 p. Zakšek, K., Podobnikar, T., Oštir, K. 2005. Solar radiation modelling. Computers & geosciences 31, 2: 233–240. Zakšek, K. 2006. Analiza vidnosti s prostorskim kotom odprtega neba. Geografski vestnik 78, 2: 97–109. Zidarič, T. 2007. Uporaba posnetkov senzorja MODIS. Diploma thesis. Ljubljana, University of Ljubljana, Faculty of civil and geodetic engineering, Departement of geodesy, Geodesy curse: 103 f. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Appendix A – Land surface temperature retrieval PHYSICAL BACKGROUND Temperature is a measure of the average kinetic energy of atomic or molecular particles in motion (Short and Stuart, 1982). The interactions between these particles lead into energy changes emitted as radiation. The thermal energy of a substance can be indicated by its kinetic (measured by a thermometer) or radiant temperature (measured by a radiometer). The Stefan- Boltzman law states the correlation between kinetic and radiant temperature: total radiation emittance (radiant temperature) integrated across all wavelengths depends on the fourth power of the absolute temperature (kinetic temperature) of a body. Therefore, the kinetic temperature (the term corresponds to the black body, otherwise brightness temperature is the appropriate term) can be measured without directly touching the body – by remote sensing. Land surface temperature (LST) can be retrieved by remote sensing that is based on the interpretation of energy measurements of a distant body. A black body is an ideal radiator on all wavelengths – it absorbs and emits the radiant energy with the maximum possible flux density, which means that all received radiation is converted into heat (Short and Stuart, 1982). According to Planck’s law (30) the measurements of the black body spectral radiance W at a certain wavelength . in the thermal infrared spectre are correlated with the kinetic temperature T of the radiating body. CC W = 1 . T = 2 2 . C . ..5 ·W . .5 ·.e.·T -1. .· ln.. 1 + ... .. C .. . 1 . -16 -2 C1 = 3.74 ·10 W m (30) C2 = 1.44 ·10-2mK In addition to temperature, the quantity of radiant emission is also controlled by the emissivity (.) of the object in the spectrum of interest. This represents the ratio of energy radiated by the material to the energy radiated by a black body at the same temperature. According to the Stefan-Boltzmann law, the radiation flux density j emitted from the Earth surface is Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. proportional (with the Stefan constant .) to the fourth power of absolute temperature T. This is true for a black body, however a real body is characterized by brightness temperature, which is a radiant temperature that is proportionally smaller than the radiant temperature of a black body for the emissivity .: 0 < . < 1 (31). j=.·. ·T 4 (31) -8 -1 -2 -4 .= 5.67 ·10Jsm K Emissivity also depends on the wavelength, however a typical engineering assumption is to assume that a surface's spectral . does not depend on it, thus . is treated as a constant (the grey body assumption). The quantity of radiation is always controlled by the kinetic temperature and . for each channel, which creates an under-determined problem for remote LST sensing because there is only one available measurement (radiance) and two parameters that need to be explained. The problem is usually solved by estimating . from other data (for example, land cover, normalised differential vegetation index – NDVI, fraction of vegetation cover – FVC, etc.). LST RETRIEVAL METHODS Satellite sensors detect the top of the atmosphere brightness temperatures, which are not the same as LST due to the atmospheric absorption and emission (figure 70). There are two “windows” in which the atmospheric influences are minimal: 3–4 µm and 8–14 µm. Approximately 80% of the energy received by the thermal sensors in the 10.5–12.5 µm wavelength region is emitted by the land surface, which makes LST relatively easy to extract from the thermal infrared signal if . is known. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. 1 4 1 3 5 2 2 Figure 70: The Sun and the Earth emit electromagnetic waves (1) that are influenced by the atmosphere (2). Part of the electromagnetic spectrum is emitted and reflected by the Earth’s surface (3). The satellite sensor records a signal that makes it possible to determine LST. The recorded LST can be used to compute AAT, which is usually measured in meteorological stations at the height of two meters (5). Slika 70: Sonce in Zemlja oddajata elektro-magnetno valovanje (EMV; 1), na katerega vpliva atmosfera (2). Na zemeljskem površju se del EMV Sonca absorbira in del odbije, nekaj EMV oddaja tudi Zemlja (3). Satelitski senzor zabeleži prejeto EMV (4) pod vplivom atmosfere (2), iz katerega lahko določimo temperaturo površja. Nato je možno iz izmerjenih vrednosti določiti še temperaturo zraka na višini dveh metrov, ki jo sicer merimo v meteoroloških opazovalnicah (5). The three prevalent methods that aim to compensate atmospheric and angular effects on LST can be used in the event of an a priori given emissivity .: single-channel methods, the split- window technique as well as the multi-angle methods (Wan and Dozier, 1996). In order to apply the last two methods at least two thermal channels are required. SINGLE CHANNEL METHOD The estimation of LST with a single thermal channel is the main advantage for single-channel methods. LST =f (p, h, z...) (32) This method requires accurate information on the vertical and horizontal distribution of air pressure p and water vapour ea in the atmosphere, solar zenith angle z, etc. Atmospheric corrections for each profile are derived from the simulated results and the known surface variables. This is a simulation of the actual conditions, thus it is the most accurate method for LST estimation. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. SPLIT-WINDOW METHOD The split-window method is based on the differential absorption in two spectral channels within a 10–12 µm atmospheric window. It is assumed that the contribution of both channels is similar and that the difference in transmittance results from the different absorption of water vapour and not due to the emissivity . differences. The split-window method is not always an optimal choice, as its coefficients a and b are strictly valid only for the data used to derive them and do not always reflect the real situation. A priori . is the critical issue for this method. A simple split-window method using channels T1 and T2 can be mathematically described with the following equation. LST = T1+ a ·(T1 - T2 )+ b (33) MULTI-ANGLE METHOD The multi-angle method is similar to the split-window method however the differential absorption occurs due to different atmospheric path-lengths when the same target is observed under different zenith angles z within the same spectral range. LST =f (.z) (34) The measurements can be made from one satellite or from two simultaneously. One of the measurements has to be made for a significantly longer path and the surface topography is also required as an input. GENERALISED SPLIT-WINDOW METHOD Another method, often used today, was developed on the basis of all described methods – it is called the generalised split-window method. . 1 -.... T1 + T2 . 1 -. ... T1 + T2 LST = C +. A + A ·+ A · .· +. B + B ·+ B · .· (35) 1232 1232 . ... 2 . ... 2 Simulations using accurate radiative transfer showed that the coefficients in the split-window method must vary with the zenith angle in order to generate a high quality LST product. The method still needs . (. is the mean . and .. is the difference in . for channels T1 and T2). The Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. coefficients (Ai, Bi, C) are obtained from the regression analysis of radiative transfer simulations (Wan and Dozier, 1996). AVAILABLE LST PRODUCTS Numerous sensors with at least one thermal infrared (TIR) channel that makes it possible to derive LST from the recorded data are currently available for Europe. The table below shows the comparison of the most interesting sensors and their LST products that are capable of recording in TIR spectre and are made available for Europe. Table 8: Comparison of the most appropriate sensors for LST measurements. Preglednica 8: Primerjava med najprimernejšimi senzorji za merjenje temperature površja. AVHRR MODIS ASTER SEVIRI Spatial resolution in nadir 1100 m V: 250 m SWIR: 500 m TIR : 1000 m V&NIR : 15 m SWIR: 30 m TIR : 90 m V: 1000 m IR: 3100 m Number of channels V: 2 SWIR: 1 TIR : 2 V: 2 SWIR: 5 TIR : 29 V&NIR : 3 SWIR: 6 TIR : 5 V&NIR: 4 SWIR: 3 TIR: 5 Swath width 3000 km 2330 km 60 km whole hemisphere (el.>10o) Satellite NOAA Terra, Aqua Terra MSG-1 Orbit 825 km high, polar, Sun synchronous 704 km high, polar, Sun synchronous 704 km high, polar, Sun synchronous geostationary (MSG-1 3.4o W) Time of equator crossing 13.55 10.30 (Terra), 13.30 (Aqua) 10.30 / Shortest time of revisit (local time) 1 day 1–2 days 16 days 15 min Where to obtain the processed LST data EOWEB (DLR) LPDAAC (NASA) LPDAAC (NASA) LandSAF A number of other sensors are appropriate for LST measurements even though there is no LST product distribution available; some interesting examples that were not investigated in detail are listed below: • AATSR on board Envisat is primarily intended for detecting the sea surface temperature; it has similar characteristics as AVHRR or MODIS, • AVHRR3 on board MetOp is primarily intended for detecting the sea surface temperature; it is an improved version of the “old” AVHRR, Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. • ETM+ on board Landsat 7 has only one IR channel with a lower temperature accuracy, but it has a 60 m spatial resolution; it is malfunctioning since 2004, and • MTI on board the MTI satellite has 11 IR channels in a 20 m spatial resolution but this is an US military satellite that is no longer operational (launched in 2000; 3 year lifespan). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Appendix B – Modelling methods MULTIPLE REGRESSION The purpose of multiple regression is to learn more about the relationship between several independent variables and a dependent variable (StatSoft, 2003). Multiple regression procedures are widely used in social and natural sciences. It allows the researcher to search for the best predictor of his problem. For example, a psychologist may wish to determine which personality variable would best predict social adjustment. The general computational problem that needs to be solved in multiple regression analysis is the fitting of a straight line to a number of points. The goal of linear regression procedures is to fit a line through the points. A line, where squared deviations of the observed points from it are minimized, has to be computed. This procedure is referred to as least squares estimation. A line in a two dimensional space can be expressed in terms of a constant c and a slope s times variable x. In the multivariate case, the regression line cannot be visualized in a two dimensional space, but can be computed just as easily. One can construct a linear equation containing all these variables; multiple regression procedures will estimate a linear equation of the following form. n y = c +.si · xi . (36) i=1 The regression line expresses the best prediction of the dependent variable y, given the independent variables xi. However, nothing is perfectly predictable and the observed points are usually scattered around the fitted regression line. The deviation of a particular point from the regression line is called the residual. The smaller the variability of the residuals around the regression line, the better the prediction. Multiple regression has certain drawbacks. The major conceptual limitation of all regression techniques is that one can only ascertain relationships, but the underlying causal mechanisms are never accurately known. Furthermore, multicollinearity is a common problem in Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. numerous correlation analyses. When a lot of variables are involved, it is often not immediately apparent that some variables are correlated – at least one of the predictor variables is redundant regarding the other variables. The fitting of higher-order polynomials of an independent variable with a mean not equal to zero can also create difficult multicollinearity problems. Most gross violations in multiple regression can be detected and should be appropriately dealt with. Outliers can seriously bias the results and this leads into biased regression coefficients. Excluding merely a single extreme case can yield a completely different set of results. CONTEXTUAL ANALYSES In GIS, raster spatial analysis is divided into four groups (Burrough and McDonnell, 1998): • local analysis only uses data from a single cell to calculate the output value, • contextual (neighbourhood) analysis performs an operation on a set of cells surrounding the selected cell within a moving window, • zonal analysis performs an operation within the zones defined by the zonal dataset, • global analysis uses all data from a raster layer to perform the operation. Contextual analysis is a part of the raster spatial analysis that creates output values for each cell location based on the value for the location and the values identified in a specified neighbourhood. They can be used to generate numerous outputs – for example: computing the relief slope, relief aspect, analytical hillshade, etc. from a digital elevation model. Even some more “simple” analyses belong to the contextual analysis – for example, “focal-mean”, “focal-median”, “focal-max”, etc. all look for a corresponding statistical measure within the “focal window”. The size and shape of the moving window depend on the analysis type: it can be a rectangle of any dimension, a circle of any radius, an annulus of any radius, and a wedge in any direction. Conceptually, a contextual analysis visits each cell in the raster and calculates the specified values with the identified neighbourhood – moving window. The value of the processed cell as well as all cell values in the identified neighbourhood are included within the computation. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The moving windows overlap, thus the cells in one moving window may also be included in the moving windows of another processed cell. WAVELET TRANSFORM Wavelet transform makes it possible to decompose a signal into several components, each of which captures the information present at a given resolution (Acerbi-Junior et al., 2006). This enables the introduction of the concept of details between various levels of resolution and if the process is inverted, the original image can be reconstructed in detail. Therefore, the wavelet transform is used today in numerous fields such as geophysics (seismic events), medicine (EKG and medical imaging), astronomy (image processing), computer science (object recognition and image compression), photography (merging images with different focus), GIS (raster generalization) and remote sensing (panchromatic sharpening). In order to understand the need for wavelets transforms one needs to understand the background of Fourier transform as well as the need for any transform at all. Most of the signals are time-domain signals in their raw format – the signals are a function of time, thus their plot is always the time-amplitude representation of the signal. However, this representation is not always the best representation of the signal, because some information can be hidden, thus the signal has to be transformed into another form in order to reveal the hidden information. Interesting information can be extracted from the frequency content of the signal: the frequency spectrum shows which frequencies can be found within the signal. The frequency content of the signal can be described by using one of the Fourier transforms – the result is usually a plot that shows how much of each frequency can be found within the signal. This is a useful way to evaluate all types of signals – for example, medical staff might not recognize a pathological heart condition from the raw ECG signal (graphical recording of heart's electrical activity), however the new computerized ECG analyzers also provide frequency information which helps them decide whether a pathological condition exists or not (Polikar, 2001). Fourier transform is a reversible transform, however only the raw or the transformed signal is available at any given time – no frequency information is available in the time-domain signal, and no time information is available in the Fourier transformed signal. This is not a problem Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. when the signal is stationary – its frequency content does not change in time, thus one does not need to know at what times a certain frequency component exists, since all frequency components exist at all times. The plot of the transformed data of a non-stationary signal, the frequency of which changes over time, can be similar – it has a few peaks, but its frequency component is not the same at all times. Therefore, the Fourier transform can be used for non- stationary signals, if one is interested merely in what spectral components can be found within the signal, but not in where they occur. However, if the time information is also required, then Fourier transform is not the appropriate transform; the time-frequency representation of the signal is needed and wavelet transform is the appropriate one to use. Another solution is short time Fourier transform, but the best it can do is an investigation of which spectral components exist at any given interval of time (Polikar, 2001). The continuous wavelet transform was developed as an alternative approach to the short time Fourier transform in order to overcome the resolution problem. Wavelet transform decomposes the signal into elementary functions – wavelets (Acerbi-Junior et al., 2006). For a one- dimensional signal f(t) the continuous wavelet transform is defined in the following equation. 1 +. . t - b . WT (a,b)= ·. f ()t ·.. .· dt (37) a -. . a . As seen in equation (37), the transformed signal is a function of two variables: scale a and translation b. The transforming function . is called the mother wavelet. Each base function is a scaled and translated version of the mother wavelet. The transformed signal actually presents the correlation between a wavelet at different scales and the signal with the scale being used as a measure of similarity. The process can be inverted and the raw signal reconstructed without any loss by (13), but it is only possible if the admissibility condition C. is defined (usually this means that the wavelet mean equals zero; Acerbi-Junior et al., 2006). 1 +.+. 1 t - b da · db f ()t = ·.. ·WT (a,b)·.. . .. · (38) C a . a . a2 .0 -. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Appendix C – Validation methods ROOT MEAN SQUARE ERROR (RMSE) RMSE is an expected value estimator of the root square of the error (Kanevski and Maignan, 2004). The error is the amount by which the estimator x differs from the quantity to be estimated x^ . The difference occurs due to the randomness or because the estimator does not account for the information that could produce a more accurate estimate. In order to define the root mean square error the mean square error (MSE) has to be defined first. MSE(x)=E((x -x^)2 ) (39) MSE is the sum of the variance and the square of the bias of the estimator. In this sense, the MSE assesses the quality of the estimator in terms of its variation and unbiasedness. The root mean squared error of n estimations is thus simply defined as the square root of the MSE. ()= RMSE x 1·.n (xi -x^)2 (40) i=j n CORRELATION (PEARSON) INDEX Correlation r is one of the most common and useful statistics. A correlation is a single number that describes the degree of relationship between two variables. This is a value between -1 and 1 that measures the degree to which two variables are linearly related (Kanevski and Maignan, 2004). If the correlation is positive, the relationship is positive (the greater the one variable, the greater the other variable); if it is negative, the relationship is negative (the greater the one variable, the smaller the other variable). n · x ·y -( x)·( y) r= ... (41) 22 22 (n ·.x -(.x))·(n ·.y -(.y)) Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. The described index is the Pearson Product Moment Correlation index; other correlations are also available – when the investigated variables are ordinal, it makes greater sense to use the Spearman Rank Correlation index, but within this thesis only the Pearson correlation index was used. AUTOCORRELATION (MORAN) INDEX The autocorrelation index is a correlation of a variable with itself, when considering observations with a gap in time (temporal autocorrelation) or space (spatial autocorrelation). If a value of an attribute on a random point is similar in the neighbouring areas, then a contiguity effect in the spatial structure exists; the phenomenon shows spatial autocorrelation (Hypergeo, 2004). Neighbouring regions tend to have identical properties or similar values with a positive autocorrelation or different qualities, i.e. alternating strong and weak values with a negative autocorrelation. When the distribution of data is random, there is no autocorrelation. The autocorrelation measures depend on the analysis scale, on the level of the grid resolution through which the distribution is observed. These measures are ratios between the covariance measured for a given interval and total variance. The most frequently used indices are the Moran index I (used within the thesis) and the Geary index. .Wij ·(xi -x)·(xj -x) I = ij (42) n ·.2 ()x TESTING THE SIGNIFICANCE OF A CORRELATION Once a correlation index is computed, one can determine the probability that the observed correlation occurred by chance by conducting a significance test. When merely a small sample is available it has to be determined that the correlation is a real one and not a chance occurrence. In this case, two exclusive hypotheses are tested: • null hypothesis: r = 0, and • alternative hypothesis: r <> 0. Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. As in all hypotheses testing, first of all the significance level has to be determined. Due to the linear correlation the degrees of freedom s equal to n–2. The significance of the correlation index is tested using the equation following a Student t-distribution (Lowry, 2007): r · n - 2 t = 1- r 2 (43) If the computed t value is greater than the t value listed in the t-distribution tables for the corresponding significance level and degrees of freedom, one can consider that the correlation index r is not random – the null hypothesis is rejected and an alternative hypothesis is accepted. CROSS-VALIDATION Cross-validation is quite different from the usual "split-sample" (calibration and validation set) method commonly used for validation (Kanevski and Maignan, 2004). In the split-sample method, only a single subset (the validation set) is used to estimate the error, instead of k different subsets. In k-fold cross-validation, one divides the data into k subsets of (approximately) equal size. The model is validated k times, each time leaving out one of the subsets from calibration, but using only the omitted subset to compute the error criterion. If k equals the sample size, this is called "leave-one-out" cross-validation. "Leave-v-out" is a more elaborate and expensive version of cross-validation that involves leaving out all possible subsets of v cases. The cross-validation is markedly superior for small data sets. SCATTERPLOT Scatterplot is a graph used in statistics to visually display and relate two quantitative variables of a multidimensional dataset by displaying the data as a collection of points, each having one coordinate on the horizontal and one on the vertical axis (Kanevski and Maignan, 2004). It is one of the basic quality control tools. A scatterplot does not require the user to specify dependent or independent variables. Either type of variable can be plotted on either axis. Scatterplots represent an association between two variables. If the pattern of dots slopes from lower left to upper right, it suggests a positive correlation between the studied variables and vice versa. A ‘best fit’ line can be drawn (the solid line within the thesis) in order to study the Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. correlation between the variables and the line “one to one”, showing the ideal relationship between the studied variables (the dashed line within the thesis). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. Appendix D – Case study areas and data sources Due to my cooperation with the Slovenian meteorological office and the German aerospace centre, AAT measurements were available only for Slovenia and Germany, thus these two areas were chosen for the case study areas. AAT interpolation and LST downscaling were performed only for Slovenia. LST to AAT parameterization was tested in both areas. Figure 71: Position of Slovenia (black) and Germany (dark grey) within Europe (data source ESRI maps). Slika 71: Položaj Slovenije (črno) in Nemčije (temno sivo) v Evropi (vir podatkov ESRI maps). Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. In this thesis numerous datasets were used with different geo-reference data, thus all data for Slovenia was projected onto a standard UTM projection (zone 33 N; WGS 84 ellipsoid); the data for Germany was projected onto the Lambert Conformal Conic projection (International 1909 ellipsoid) with the following parameters. Table 9: Projection parameters for Germany. Preglednica 9: Parametri projekcije za Nemčijo. 1st standard parallel 48°40’ 2nd standard parallel 53°40’ Central meridian 10°30’ Latitude of projection’s origin 51° False easting 0 m False northing 0 m Due to the problems with cloudy weather, only some dates were appropriate for testing the parameterization. The following dates were used for testing in Slovenia (at 6:00, 13:00 or 21:00 GMT – not more than 3 examples per day): 30.05.2005 28.07.2005 27.10.2005 03.06.2005 01.09.2005 28.10.2005 08.06.2005 02.09.2005 29.10.2005 19.06.2005 13.09.2005 09.11.2005 20.06.2005 14.09.2005 10.11.2005 23.06.2005 15.09.2005 18.11.2005 03.07.2005 26.09.2005 18.12.2005 26.07.2005 28.09.2005 21.12.2005 27.07.2005 30.09.2005 24.12.2005 The following dates were used for testing in Germany (every hour between 4:00 and 21:00 GMT – not more than 18 examples per day): 03.06.2005 28.07.2005 22.09.2005 21.06.2005 24.08.2005 13.10.2005 23.06.2005 29.08.2005 28.10.2005 19.07.2005 31.08.2005 30.10.2005 21.07.2005 05.09.2005 08.11.2005 27.07.2005 08.09.2005 12.11.2005 Zakšek, K. 2007. Air temperature determination from land surface temperature model of high spatial and temporal resolution. PhD thesis. Ljubljana, UL, FGG, Department of Geodetic Engineering, Postgraduate study. New algorithms were developed within the study by using data stored in various formats, thus the following software was used: • RSI IDL 6.3 (to model and to plot diagrams), • GNU Octave 2.9 (to model solar radiance exposure), • ESRI ArcGIS 9.1 (to prepare relief and land cover data and to map the results), • MS Access 2003 (to organize the meteorological measurements into a database), • MS Excel 2003 (to plot diagrams), • MODIS Reprojection Tool 3.3 (to prepare MODIS data), • NCSA HDFView 2.3 (to visualize data stored in HDF files), and • BEAM 3.6. (to test AVHRR and AATSR data). Data provided by the following institutions was used in the modelling procedures: • © Deutscher Wetterdienst, 2005–2007, • © EEA, 2005, • © Esri, 2004, • LandSAF, • NASA, • © Surveying and Mapping Authority of the Republic of Slovenia, 2005–2007, • Environmental Agency of the Republic of Slovenia.