Interpolacja przestrzenna należy do grupy metod analizy przestrzen-nej opartych na transformacji danych (Strobl, 1999; Urbański, 2008). Umoż-liwia ona oszacowanie wielkości analizowanej zmiennej w dowolnym punk-cie obszaru badań na podstawie dyskretnej, najczęśpunk-ciej punktowej informacji, pozyskanej poprzez pomiar lub obserwację. Procedurę taką można także na-zwać modelowaniem pola zmiennej, choć model, opisując strukturę pola, nie zawsze musi prowadzić do wytworzenia informacji ciągłej przestrzennie (Szymanowski, Kryza, 2011b).
Pojęcie interpolacji jest w analizie przestrzennej używane w szer-szym sensie, niż gdyby wywodzić je ze ściśle matematycznej definicji, mówiącej o jednoznacznym wyznaczaniu funkcji, przechodzącej w danym przedziale (tu: obszarze) przez znane punkty węzłowe. Po pierwsze, metody interpolacji przestrzennej obejmują także algorytmy przybliżające (wygładza-jące), a więc takie, które nie wymuszają przeprowadzenia modelowanej po-wierzchni ściśle przez punkty obserwacji. W takim ujęciu interpolacja jest szczególnym przypadkiem aproksymacji funkcyjnej, której naczelną zasadą jest minimalizacja powstałych błędów aproksymacji (Szymanowski i in., 2006). Po drugie, w sytuacjach, gdy opróbowanie nie obejmuje całego anali-zowanego regionu, szacowanie wielkości zmiennej odnosi się do punktów leżących poza zakresem przestrzennym danych wejściowych, którego przy-bliżeniem może być otoczka wypukła danych (convex hull), a więc jest de facto przestrzenną ekstrapolacją (Burrough, McDonnell, 1998). Zatem podstawowy cel interpolacji przestrzennej może być zdefiniowany jako
sza-cowanie (estymacja) wartości cechy w dowolnej lokalizacji obszaru badań, także poza zakresem przestrzennym danych obserwacyjnych. W literaturze przedmiotu termin estymacja jest używany jako synonim przestrzennej inter-polacji, zwłaszcza w kontekście stosowania metod geostatystycznych. Poję-cie interpolacji przestrzennej, jako metody obejmującej szacowanie (estyma-cję) czy też przybliżanie (aproksyma(estyma-cję) wartości zmiennej, zarówno we-wnątrz (interpolacja), jak i na zewe-wnątrz (ekstrapolacja) obszaru badań jest powszechnie akceptowane w analizie przestrzennej i GIS (Lam, 1983; Mitas, Mitasova, 1999; Dubois, Galmarini, 2004). Czasem bywa także zastępowane terminami: spatial prediction – prognoza przestrzenna (Hengl, 2007) lub spa-tialization (Tveito, 2007; COST Action 719 Final Report, 2008). Ten ostatni, w wersji spolszczonej, jako „spacjalizacja”, od niedawna pojawia się w pol-skiej literaturze przedmiotu (Ustrnul, Czekierda, 2003).
Niezależnie od przyjętej nazwy, głównym celem interpolacji jest otrzymanie ciągłej przestrzennie informacji o wartości zmiennej oraz zapew-nienie możliwości jej graficznej prezentacji w formie mapy. Zagadzapew-nienie interpolacji w kartografii ma bardzo długą historię i jest nierozerwalnie zwią-zane z opracowywaniem map izoliniowych, które do niedawna były dominu-jącą formą prezentacji ciągłych zjawisk klimatologicznych na mapach.
W metodzie tej problematycznym pozostaje kształt izolinii i sposób jej prze-prowadzania pomiędzy punktami wyznaczanymi na osiach interpolacyjnych, gdyż aparat matematyczny nie uwzględnia rzeczywistych cech środowisko-wych interpolowanej zmiennej. Częściowym rozwiązaniem tego problemu była koncepcja tzw. interpolacji geograficznej z autorską modyfikacją prze-biegu izolinii i dopasowywaniu ich do charakteru rozkładu zjawiska (Gor-czyński, 1918; Romer, 1938, 1948/1949). Proces taki, pomimo iż może uwzględniać dodatkowe informacje pominięte w danych wejściowych, a tym samym podnosić wartość mapy, jest obarczony wysoką dozą subiektywizmu.
Ceną za bardziej „naturalny” przebieg izolinii jest obniżenie matematycznej precyzji interpolacji (Mościbroda 1999).
Współczesne techniki aproksymacji danych przestrzennych dążą do obiektywizmu i powtarzalności procedur, które pozwalają na ocenę jako-ści (niepewnojako-ści) dokonanej interpolacji oraz zapewniają możliwość doboru optymalnej metody w zależności od charakteru zjawiska i przestrzennej dys-trybucji informacji źródłowej (Szymanowski i in., 2006). W wielu z tych technik wykorzystywana jest dodatkowa informacja środowiskowa, podno-sząca precyzję i jakość interpolacji (interpolacja wielowymiarowa). Należy zaznaczyć, iż mimo dążenia do obiektywizmu i wykorzystania ilościowego, matematycznego lub statystycznego aparatu teoretycznego, interpolacja prze-strzenna nie jest jedynie semi-automatyczną procedurą wykonywaną w śro-dowisku cyfrowym komputera. Decyzje podejmowane przez wykonawcę wymagają wiedzy eksperckiej, zarówno o naturze wykorzystywanych metod, jak i o właściwościach i zależnościach przestrzennych interpolowanej zmien-nej. Podsumowując, obecnie tworzenie modeli/map ma często charakter ana-lityczny i polega na obliczaniu wielkości zmiennej w każdym punkcie anali-zowanego obszaru za pomocą programów komputerowych, na podstawie pomiarów lub obserwacji terenowych oraz uzupełniającej informacji środo-wiskowej (Hengl, 2007).
Metody interpolacji przestrzennej zostały wypracowane dla zastoso-wań w różnych dyscyplinach i z przeznaczeniem dla różnorodnych zbiorów danych i zmiennych. Takie czynniki, jak: wielkość zbioru wejściowego, roz-mieszczenie przestrzenne danych, cel interpolacji i charakter interpolowanej zmiennej wpływają na wybór odpowiedniej metody. Zagadnienie znalezienia optymalnej metody dla danego zbioru danych przestrzennych jest kluczowe, ale często trudne do jednoznacznego rozwiązania. Pomocne są tu różnorodne klasyfikacje, grupujące metody w zależności od wybranych cech i założeń
teoretycznych. Należy jednak pamiętać, iż szereg metod może dawać bardzo podobne rezultaty, mimo że ich fundamenty teoretyczne są krańcowo różne (COST Action 719 Final Report, 2008).
Za idealną można uznać sytuację, gdy rozkład przestrzenny analizo-wanej zmiennej jest determinowany przez skończony zbiór czynników, pod-legających znanym prawom fizyki. W takim przypadku formuła interpolacji jest znana, a prognoza dla każdego punktu obszaru odbywa się jednoznacz-nie, za pomocą modelu deterministycznego (Hengl, 2007). Jednak w rzeczy-wistości relacje pomiędzy zmienną a środowiskiem są tak skomplikowane, iż nie może być ona modelowana wprost (Heuvelink, Webster, 2001), a mo-del przyjmuje postać częściowo lub całkowicie niedeterministyczną (Hengl, 2007). Ten fakt jest uwzględniony w podstawowym podziale metod interpo-lacyjnych na deterministyczne i probabilistyczne oraz ich kombinacje (Ustrnul, Czekierda, 2003; COST Action 719 Final Report, 2008). Do grupy metod deterministycznych zalicza się m.in. takie metody, jak: metoda od-wrotnej ważonej odległości (inverse distance weighting – IDW) czy funkcje wielomianowe (funkcje sklejane, radialne funkcje bazowe, trend powierzch-niowy). Do opartych na koncepcji losowości metod geostatystycznych zali-czane są podstawowe techniki krigingowe: kriging zwyczajny (ordinary kri-ging – OK), krikri-ging prosty i kokrikri-ging. Przy takim podziale pewien problem z zaklasyfikowaniem stwarzają metody regresji, jak liniowa regresja wielo-czynnikowa (multiple linear regression – MLR). Model regresji jest zalicza-ny do metod stochastyczzalicza-nych przy założeniu braku autokorelacji przestrzen-nej i normalnym rozkładzie reszt regresji. Jeśli jednak jedynym założeniem modelu jest możliwość interpretacji procesu przestrzennego prawami fizycz-nymi, to można go zaklasyfikować jako technikę deterministyczną (COST Action 719 Final Report, 2008). Do ostatniej grupy, metod złożonych (kom-binowanych), a więc takich, które składają się z komponentów
determini-stycznego i stochadetermini-stycznego zalicza się kriging uniwersalny, kriging z ze-wnętrznym dryftem (przesunięciem) i kriging resztowy (residual kriging, regression kriging – RK).
Pozostałe podziały algorytmów interpolacyjnych pozwalają na wy-różnienie metod:
globalnych i lokalnych, w zależności od tego, czy budowa modelu oparta jest na całej populacji danych obserwowanych, czy też na jej podzbiorach,
wiernych i wygładzających, z uwagi na możliwość wiernego zachowania danych pomiarowych lub tendencję do wygładzania estymowanej po-wierzchni,
jedno- i wielowymiarowych, w zależności od tego, czy w procedurze in-terpolacji wykorzystywane są dodatkowe zmienne objaśniające, zwane także, z uwagi na powiązania przyczynowo-skutkowe i naturalne uwarun-kowania, zmiennymi środowiskowymi.
Obecnie dla wykonania przestrzennej interpolacji wykorzystuje się co najmniej kilkadziesiąt różnych algorytmów. Przykładowo, w opracowaniu Li i Heapa (2008) omówiono aplikacje 62 metod i ich odmian, zamieszczone w 51 opracowaniach dotyczących interpolacji różnych elementów środowi-ska przyrodniczego. Powstają pytania:
które z tych metod można uznać za odpowiednie w przypadku temperatu-ry powietrza?
czy można ich używać dla każdego obszaru badań, skali przestrzennej i stopnia agregacji danych, niezależnie od liczby i rozmieszczenia punk-tów danych?
czy istnieje metoda uniwersalna, którą można stosować do estymacji tem-peratury powietrza?
Poszukiwanie odpowiedzi na takie pytania przewija się w literaturze przedmiotu, zwłaszcza w okresie ostatnich 20 lat. Porównanie wyników ze wszystkich dostępnych metod jest praktycznie niemożliwe, a sytuację komplikuje szereg czynników związanych z celem opracowania, charakterem zbioru danych, pożądaną dokładnością i rozdzielczością wyników, a także z naturą zjawiska w konkretnym obszarze badań. Większość tych zagadnień znalazło podsumowanie w ramach europejskiego programu badawczego COST 719 „The Use of GIS in Climatology and Meteorology”, zwłaszcza w pracach grupy roboczej WG 2 „Spatial Interpolation”. Aktywność tej gru-py, w kontekście spacjalizacji danych meteorologicznych i klimatologicz-nych, koncentrowała się na inwentaryzacji metod i narzędzi interpolacyjnych zaimplementowanych w środowisku GIS i oprogramowaniu statystycznym.
Projekt COST 719 był realizowany w latach 2001–2006 i został podsumowa-ny na konferencji w Grenoble w 2006 r. i w dwóch wydawnictwach książ-kowych (Dobesch i in., 2007; COST Action 719 Final Report, 2008). Dodat-kowo zostały zorganizowane dwie konferencje pod tytułem „Spatial Interpo-lation in Meteorology and Climatology” w Budapeszcie w 2004 (Szalai i in., 2007) i 2009 r., na których zostały przedstawione także prace naukowców niezwiązanych z projektem COST 719.
Spektrum metod używanych do interpolacji temperatury powietrza jest bardzo szerokie, co znajduje odzwierciedlenie w literaturze przedmiotu, obejmującej wiele opracowań różniących się skalą i rozdzielczością prze-strzenną, położeniem i cechami obszaru badań, stopniem agregacji danych (rozdzielczość czasowa), liczbą i rozmieszczeniem stacji meteorologicznych.
Wielkość obszarów analizy jest bardzo różna i zmienia się od skali globalnej (Hijmans i in., 2005) i kontynentalnej (Willmott, Matsuura, 1995; Hong i in., 2005; Daly i in., 2008; Hutchinson i in., 2009; Hogewind, Bissolli, 2011;
Krahenmann i in., 2011) do skali lokalnej, gdzie powierzchnia obszaru
inter-polacji to zaledwie 300–500 km2 (Daly i in., 2007; Szymanowski, 2004;
Szymanowski, Kryza, 2009a, b, 2012a). Z wielkością obszaru analizy zwią-zana jest często także rozdzielczość przestrzenna – z reguły duża dla obsza-rów o niewielkiej powierzchni, np.: 25 m (Szymanowski, Kryza, 2009, 2012a) czy 2” (Daly i in., 2007), i mała dla rozległych połaci kontynentów, np. 300” (Hutchinson i in., 2009; Hogewind, Bissolli, 2011).
Liczba opisanych w literaturze metod, które były wybierane do te-stów i porównania wyników interpolacji temperatury powietrza sięga kilku-dziesięciu i obejmuje znane metody podstawowe (np. IDW, OK, MLR, RK), jak i takie, które mają szczególne znaczenie dla estymacji elementów klima-tu, na przykład tzw. metoda gradientowa (Lennon, Turner, 1995; Willmott, Matsura, 1995; Courault, Monestiez, 1999; Stahl i in., 2006), metoda TPS (thin-plate splines) z grupy funkcji sklejanych (Hutchinson, 1989, 1995; Hu-tchinson, Bischof, 1983) lub metody: AURELHY (Benichou, 1986; Beni-chou, Le Breton, 1987), PRISM (Daly i in., 1994, 2002, 2008) i MISH (Szen-timrey i in., 2005; Szen(Szen-timrey, Bihari, 2007; Szen(Szen-timrey i in., 2007, 2011), które zostały rozwinięte celowo na potrzeby spacjalizacji klimatologicznej.
Temperatura powietrza, a zwłaszcza sterujące nią w przestrzeni me-chanizmy fizyczne i uwarunkowania środowiskowe są stosunkowo dobrze rozpoznane. Wydaje się więc, że niezależnie od tego czy algorytm interpola-cyjny oparty jest na koncepcjach matematycznych czy statystycznych, głów-nym kryterium zasadności doboru metody jest jej jedno- lub wielowymiaro-wość. Już do samych testów opisanych w publikacjach metody wielowymia-rowe są wybierane znacznie częściej niż jednowymiawielowymia-rowe. Można także uznać za prawidłowość, że w sytuacji przeprowadzenia porównania wyników uzyskanych metodami jedno- lub wielowymiarowymi, ostatecznie wybrana zostanie jedna z tych ostatnich. Przeciwne przypadki są niezmiernie rzadkie (np. Chuanyan i in., 2005; Serbin, Kucharik, 2009). Także w opracowaniach,
w których nie przeprowadzano porównania metod, a obligatoryjnie wskazy-wano optymalną metodę, zawsze decydowskazy-wano się na technikę wielowymia-rową, np. liniową lub nieliniową metodę regresji (White, 1979; Gomez i in., 2008; Esteban i in., 2009; Guan i in., 2009; Hiebl i in., 2009; Boi i in., 2011), kriging resztowy (Courault, Monestiez, 1999; Tveito i in., 2000; Bjornsson i in., 2007; Perčec Tadić, 2010; Tietavainen i in., 2010) lub TPS (Hijmans i in., 2005; Hong i in., 2005; Hutchinson i in., 2009). Podsumowując, można stwierdzić, iż wśród metod uznawanych w literaturze za optymalne, dominu-ją te, które uwzględniadominu-ją w estymacji rolę wysokości bezwzględnej oraz in-nych predyktorów (m.in. MLR, TPS, RK, kombinacja MLR+IDW). Stosun-kowo mało jest także prac, w których interpolację wykonano metodami roz-winiętymi i dedykowanymi do zastosowań w klimatologii (np. PRISM, MISH). Z jednej strony są to techniki stosunkowo nowe (MISH), które wy-magają dostępu do specjalnego oprogramowania (MISH, ANUSPLIN) i za-awansowanej wiedzy do parametryzacji metody (PRISM). Z drugiej strony w literaturze nie przedstawiono przekonujących dowodów na ich przewagę nad bardziej znanymi i rozpowszechnionymi metodami, zaimplementowa-nymi we wiodących pakietach oprogramowania statystycznego i GIS.
Powyższe prawidłowości potwierdzają także opracowania wykonane w ramach akcji COST 719. W zdecydowanej większości aplikacji dla danych termicznych, uśrednionych na poziomie lat i miesięcy decydowano się na zastosowanie metod wielowymiarowych. Jako potencjalne predyktory rozkładu temperatury wykorzystywano głównie wysokość bezwzględną, współrzędne geograficzne oraz odległość od morza (Szalai i in., 2007; Tvei-to, 2007; Ustrnul, Czekierda, 2003). Najczęściej stosowaną metodą był kri-ging resztowy, za pomocą którego wykonano mapy średnich rocznych i mie-sięcznych wartości temperatury dla Hiszpanii, Fennoskandii, Słowenii i Pol-ski (COST Action 719 Final Report, 2008; Ustrnul, Czekierda, 2003, 2005).
W Portugalii, oprócz RK, dla niektórych miesięcy, stosowano kokriging (COST Action 719 Final Report, 2008), a dla Niemiec wybrano metodę re-gresji wieloczynnikowej (Klimaatlas Bundesrepublik Deutschland, 2001).
Na Węgrzech z kolei, preferowano opracowaną tam, dedykowaną danym klimatologicznym metodę MISH (Szalai i in., 2007).