UNIVERZA V LJUBLJANI FILOZOFSKA FAKULTETA ODDELEK ZA GEOGRAFIJO
MARUŠA ŠKERL
Uporaba satelitskih posnetkov Landsat in Sentinel za preučevanje širjenja pozidanih površin v Mestni občini
Ljubljana
Diplomsko delo
Ljubljana, 2021
UNIVERZA V LJUBLJANI FILOZOFSKA FAKULTETA ODDELEK ZA GEOGRAFIJO
MARUŠA ŠKERL
Uporaba satelitskih posnetkov Landsat in Sentinel za preučevanje širjenja pozidanih površin v Mestni občini
Ljubljana
Diplomsko delo
Mentor: izr. prof. dr. Marko Krevs Študijski program: Geografija - E
Ljubljana, 2021
Zahvala
Rada bi se zahvalila mentorju izr. prof. dr. Marku Krevs za vse predloge, ideje in spodbude
med pisanjem, hkrati pa tudi asist. Alešu Grlju za ves čas in pomoč pri praktičnem delu naloge. Navsezadnje pa bi se zahvalila tudi moji družini, ki je verjela vame in mi stala ob
strani ob napornih dneh.
Izvleček
Uporaba satelitskih posnetkov Landsat in Sentinel za preučevanje širjenja pozidanih površin v Mestni občini Ljubljana
V diplomski nalogi smo zaznavali pozidane površine v Mestni občini Ljubljana s tremi različnimi sateliti med leti 2000, 2013 in 2020. Uporabili smo tri različne metode: nadzorovano klasifikacijo, spektralne indekse z mejo po Otsujevi metodi in spektralne indekse z ročno nastavljeno mejo. Namen naloge je bil torej ugotoviti, katera izmed treh preprostih metod daljinskega zaznavanja daje najboljše rezultate. Ugotovili smo, da je to nadzorovana klasifikacije, prav tako pa sta tudi Normalized Difference Vegetation Index (NDVI) in Built- up Index (BU) z ročno nastavljeno mejo prikazala dokaj dober rezultat. Najslabše rezultate je dala metoda spektralnih indeksov z mejo po Ostujevi metodi. Hkrati smo potrdili izkušnje avtorjev, ki so se ukvarjali s podobnimi tematikami, da težko primerjamo satelitske slike večih satelitov, zaradi različne kakovosti (Sabo, Corbane, Ferri, 2017).
Ključne besede: pozidane površine, daljinsko zaznavanje, Mestna občina Ljubljana, satelitski posnetki, klasifikacija
Abstract
Use of Landsat and Sentinel satellite images for studying the spread of built-up areas in the Municipality of Ljubljana
In this paper we detected built-up areas in the Municipality of Ljubljana with three different satellites between years 2000, 2013 and 2020. We used three different methods: supervised classification, spectral indices with a limit according to the Otsu method and spectral indices with a manually set limit. The purpose of the task was to determine which method gives the best results. We found that it is a supervised classification, and also Normalized Difference Vegetation Index and Built-up Index with manually set limit which showed a fairly good result.
The worst results were obtained by the method of spectral indices with a limit according to the Ostuje method. At the same time, we agree with other authors which have dealt with similar topics that it is difficult to compare satellite images of multiple satellites, due to different quality (Sabo, Corbane, Ferri, 2017).
Key words: Built-up areas, remote sensing, Municipality of Ljubljana, satellite images, classification
1
Kazalo
1. Uvod ... 4
1. 1. Namen, cilji, hipoteza ... 4
2. Teoretična izhodišča ... 5
2. 1. Ključni pojmi ... 5
2. 2. Pregled literature ... 5
2. 3. Preučevano območje ... 6
3. Metodologija ... 7
3. 1. Priprava grafičnih podatkov ... 7
3. 1. 1. Vektorski podatki ... 8
3. 1. 2. Satelitske slike ... 8
3. 2. Izbrani spektralni indeksi ... 11
3. 2. 1. NDVI indeks ... 11
3. 2. 2. NDBI indeks ... 12
3. 2. 3. UI indeks ... 12
3. 2. 4. BU indeks ... 13
3. 3. Ročna metoda nastavljanja meje ... 13
3. 4. Otsujeva metoda nastavljanja meje ... 14
3. 5. Nadzorovana klasifikacija ... 14
3. 6. Ocena kakovosti klasifikacije ... 15
3. 7. Primerjava posnetkov po časovnih obdobjih ... 18
3. 8. Uporabljena orodja ... 19
4. Rezultati ... 20
4. 1. Nadzorovana klasifikacija ... 20
4. 2. Spektralni indeksi z mejo po Otsujevi metodi ... 22
4. 2. 1. NDVI ... 22
4. 2. 2. NDBI ... 23
4. 2. 3. UI ... 24
4. 2. 4. BU ... 25
4. 3. Spektralni indeksi z ročno nastavljeno mejo ... 26
4. 3. 1. NDVI ... 28
4. 3. 2. NDBI ... 29
4. 3. 3. UI ... 30
4. 3. 4. BU ... 31
5. Sklep ... 33
5. Summary ... 35
2
6. Viri in literatura ... 37
Kazalo slik
Slika 1: Preučevano območje. ... 6Slika 2: Izbrane metode zaznavanja pozidanih površin. ... 7
Slika 3: Primer vstavljene formule za indeks NDVI v orodju Raster calculator. ... 11
Slika 4: Uporaba Swipe Layer orodja. ... 14
Slika 5: Primer učnih vzorcev za leto 2000. ... 15
Slika 6: Točke, ki smo jih uporabili za oceno kakovosti klasifikacije za leto 2000. ... 16
Slika 7: Nadzorovana klasifikacija pred reklasifikacijo. ... 20
Slika 8: Pozidane površine v Mestni občini Ljubljana izračunane z metodo nadzorovane klasifikacije. ... 21
Slika 9: Pozidane površine v Mestni občini Ljubljana izračunane z NDVI indeksom in mejo nastavljeno po Otsujevi metodi. ... 23
Slika 10: Pozidane površine v Mestni občini Ljubljana izračunane z NDBI indeksom in mejo nastavljeno po Otsujevi metodi. ... 24
Slika 11: Pozidane površine v Mestni občini Ljubljana izračunane z UI indeksom in mejo nastavljeno po Otsujevi metodi. ... 25
Slika 12: Pozidane površine v Mestni občini Ljubljana izračunane z BU indeksom in mejo nastavljeno po Otsujevi metodi. ... 26
Slika 13: Pozidane površine v Mestni občini Ljubljana izračunane z NDVI indeksom in ročno nastavljeno mejo. ... 28
Slika 14: Pozidane površine v Mestni občini Ljubljana izračunane z NDBI indeksom in ročno nastavljeno mejo. ... 29
Slika 15: Pozidane površine v Mestni občini Ljubljana izračunane z UI indeksom in ročno nastavljeno mejo. ... 30
Slika 16: Pozidane površine v Mestni občini Ljubljana izračunane z BU indeksom in ročno nastavljeno mejo. ... 31
Kazalo preglednic
Preglednica 1: Podatki o izbranih satelitskih slikah. ... 8Preglednica 2: Značilnosti satelita Landsat 5. ... 9
Preglednica 3: Značilnosti satelita Landsat 8. ... 9
Preglednica 4: Značilnosti satelita Sentinel 2. ... 10
Preglednica 5: Formula NDVI indeksa. ... 12
Preglednica 6: Formula NDBI indeksa. ... 12
Preglednica 7: Formula UI indeksa. ... 12
Preglednica 8: Formula BU indeksa. ... 13
Preglednica 9: Interpretacija kappa koeficienta. ... 15
Preglednica 10: Uvoz atributivnih tabel v Microsoft Exel. ... 17
Preglednica 11: Izračun vrednosti p_e. ... 17
Preglednica 12: Izračun vrednosti p_o. ... 18
Preglednica 13: Izračun kappa indeksa. ... 18
3
Preglednica 14: Delež pozidanih površin in kappa koeficient za nadzorovano klasifikacijo. . 22 Preglednica 15: Delež pozidanih površin in kappa koeficient za indeks NDVI in mejo nastavljeno po Otsujevi metodi. ... 23 Preglednica 16: Delež pozidanih površin in kappa koeficient za indeks NDBI in mejo nastavljeno po Otsujevi metodi. ... 24 Preglednica 17: Delež pozidanih površin in kappa koeficient za indeks UI in mejo nastavljeno po Otsujevi metodi. ... 25 Preglednica 18: Delež pozidanih površin in kappa koeficient za indeks BU in mejo nastavljeno po Otsujevi metodi. ... 26 Preglednica 19: Delež pozidanih površin in kappa koeficient za indeks NDVI in ročno nastavljeno mejo. ... 28 Preglednica 20: Delež pozidanih površin in kappa koeficient za indeks NDBI in ročno nastavljeno mejo. ... 29 Preglednica 21: Delež pozidanih površin in kappa koeficient za indeks UI in ročno nastavljeno mejo. ... 30 Preglednica 22: Delež pozidanih površin in kappa koeficient za indeks BU in ročno nastavljeno mejo. ... 32
Kazalo grafov
Graf 1: Ročno nastavljene mejne vrednosti med pozidanimi in nepozidanimi površinami na klasifikacijah s spektralnimi indeksi. ... 27
4
1. Uvod
Daljinsko zaznavanje je ena izmed uporabnih metod pridobivanja informacij o Zemljinem površju. Preučevanje pozidanih površin v mestih je ena izmed pogosto zastopanih kategorij zaznavanja površja, saj se dotika človeka in njegovega delovanja v prostoru. Ker so mesta večje homogene enote pozidanih površin, nam njihovo zaznavanje omogočajo satelitski posnetki srednjih prostorskih ločljivosti, ki so široko dostopni na svetovnem spletu.
1. 1. Namen, cilji, hipoteza
Namen diplomske naloge je s pomočjo satelitskih posnetkov Landsat in Sentinel ter preprostih klasifikacijskih postopkov prikazati širjenje pozidanih površin med leti 2000, 2013 in 2020 v Mestni občini Ljubljana in primerjati uspešnost različnih metod.
Cilji:
• Poiskati primerne satelitske slike za tri časovne prereze, in sicer za leto 2000, 2013 in 2020.
• Izbrati primerne spektralne indekse za pozidane površine in primerno nadzorovano klasifikacijo.
• Spektralnim indeksom ustrezno nastaviti mejo med pozidanimi in nepozidanimi površinami na dva načina – ročno in z Otsujevo metodo.
• Preveriti kakovost klasifikacij s kappa koeficientom.
• V obliki kart prikazati dobljene rezultate in jih interpretirati.
Hipoteze:
• Klasifikacije bodo konkretno nakazovale, da so se pozidane površine v Mestni občini Ljubljana skozi celotno časovno obdobje povečale.
• Metoda zaznavanja pozidanih površin z spektralnimi indeksi bo nasploh dala boljše rezultate kot metoda nadzorovane klasifikacije.
• Built-up indeks (BU) bo, zaradi večje kompleksnosti, podal boljše rezultate kot ostali spektralni indeksi.
• Ocena kakovosti klasifikacije po kappa koeficientu bo vsaj pri eni klasifikaciji nad 0.8.
5
2. Teoretična izhodišča
2. 1. Ključni pojmi
Za splošno razumevanje naloge moramo na začetku razložiti nekaj ključnih pojmov, ki se pojavljajo skozi celotno nalogo.
Daljinsko zaznavanje je zaznavanje odbitega elektromagnetnega valovanja od površja, kar lahko uporabimo v analizah geografskih informacijskih sistemov. Zaznavamo lahko atmosferske pojave, vodne pojave in kopenske pojave na Zemlji, nekatere študije pa so usmerjene tudi v preučevanje površja drugih planetov (Oštir, 2006). Glavni metodi zaznavanja pozidanih površin v diplomski nalogi sta uporaba spektralnih indeksov in nadzorovana klasifikacija. Pri nadzorovani klasifikaciji gre za metodo, kjer na satelitskem posnetku sami izrišemo učne vzorce in s tem program naučimo, kakšne vrednosti ima določen tip rabe tal (Oštir, 2006). Spektralni indeksi so preprosti matematični izračuni z izbranimi spektralnimi pasovi posameznega satelita (Zha, Gao, Ni, 2003). Gre za enega izmed postopkov, ki dane informacije pretvori v nove s pomočjo spektralnih pasov. Spektralni pas je »merilo, ki določa razpon valovnih dolžin v opazovanem signalu« (Oštir, 2006, str. 222). Več kot ima senzor spektralnih pasov, boljša je spektralna ločljivost. Poleg spektralne ločljivosti poznamo še prostorsko, radiometrično in časovno. Rezultat izvedenih metod je torej klasifikacija, kjer gre za »kvantitativni sloj, ki ga kasneje lahko uporabimo v analizah s sistemi GIS« (Oštir, 2006, str. 175).
2. 2. Pregled literature
Od domače literature smo se pri teoretičnih osnovah največkrat opirali na Oštirjevo delo z naslovom Daljinsko zaznavanje (2006), pri praktičnem delu pa na predhodne diplomske in magistrske naloge s podobno tematiko iz Fakultete za gradbeništvo in geodezijo ter Filozofske fakultete, Oddelek za geografijo.
Ob pregledu tuje literature smo ugotovili, da je bilo tovrstnih raziskav o daljinskem zaznavanju pozidanih površin svetovno izvedenih že veliko. Najverjetneje zato, ker gre za razmeroma preproste metode, ki lahko podajo kakovostne rezultate o Zemljinem površju. Predvsem se raziskave oklepajo spektralnih indeksov, tako kot mi, s to razliko, da se pojavljajo tudi drugi spektralni indeksi, ki niso bili uporabljeni v tej nalog. Tako avtorji Setiyonoki, Hastuti in Widayani navajajo tudi spektralna indeksa na vidni osnovi, in sicer VrNIR-BI in VgNIR-BI (2017, str. 2). Tuja literatura se v tovrstnem daljinskem zaznavanju v večini navezuje na satelitske posnetke Landsat in Sentinel, kar potrjuje primernost izbranih posnetkov za to nalogo.
Poleg zaznavanja pozidanih površin veliko člankov kot osrednjo temo predstavlja zaznavanje vodnih površin in s tem povezanih poplav, gozdnih površin zaradi krčenja tropskega deževnega gozd in drugo.
Veliko slabost domače kot tudi tuje literature smo prepoznali v pomanjkljivem zapisu konkretnega dela v izbranih programih, torej opisi uporabe različnih orodij, zato smo to nalogo prilagodili temu, da bi sami kar se da natančno opisali postopke dela.
6 2. 3. Preučevano območje
Izbrano območje preučevanja je Mestna občina Ljubljana, ki leži v osrednji Sloveniji in je občina z največ prebivalci v državi. V zadnjih nekaj letih močno narašča proces suburbanizacije in s tem tudi gradnja novih stanovanjskih naselij in vse potrebne infrastrukture na obrobju mesta, zato je pričakovati trend širjenja pozidanih površin. Ker gre za osrednjo občino Slovenije, smo mnenja, da bodo razlike v letih tu najbolj izrazite.
Teoretično v nalogi smatramo za pozidane površine vse površine, ki so prekrite z gradbenim materialom, torej vse stavbe, ceste, železniški tiri in parkirišča.
Slika 1: Preučevano območje.
7
Izbrane metode zaznavanja pozidanih površin
3. Metodologija
Celostna metodologija zajema več korakov dela. Od priprave vseh vektorskih podatkov in satelitskih slik, izbire spektralnih indeksov in nadzorovane klasifikacije, postavljanja meje med pozidanimi in nepozidanimi površinami, do ocene klasifikacij in interpretacije rezultatov. Delo smo opravljali v različnih orodjih, kar je natančneje predstavljeno v poglavju 3. 8.
Uporabili smo dva pristopa zaznavanja pozidanih površin. Prvi se je nanašal na spektralne indekse za pozidane površine. Tem smo v nadaljevanju nastavili mejo med pozidanimi in nepozidanimi površinami z Otsujevo metodo in ročno. Drugi pristop je bila nadzorovana klasifikacija, pri kateri smo vnaprej določili območja, za katera želimo, da jih program zazna kot pozidana. V diplomski nalogi smo torej pozidane površine zaznavali s tremi različnimi metodami.
Slika 2: Izbrane metode zaznavanja pozidanih površin.
3. 1. Priprava grafičnih podatkov
Za doseganje vseh zadanih ciljev smo zbrali ustrezne podatke: satelitske slike in vektorske podatke meje Mestne občine Ljubljana. V nadaljevanju smo morali uskladiti tudi koordinatne sisteme posameznih slojev. Tako smo v programu Arc Map 10. 8. 1. vektorskemu sloju z ukazom Project spremenili koordinatni sistem iz D48_Slovenia_TM v WGS_1984_UTM_Zone_33N, da se je ujemal s satelitskimi slikami.
Spektralni indeksi za pozidane površine
Nadzorovana klasifikacija
Ročno nastavljena
meja
Meja nastavljena
z Otsujevo metodo
Metoda največjega
verjtja
8 3. 1. 1. Vektorski podatki
Izbrani vektorski podatki, ki smo jih potrebovali so meje Mestne občine Ljubljana. Podatke sem pridobili na spletni strani Geodetske uprave Republike Slovenije.
3. 1. 2. Satelitske slike
Ena izmed pomembnejših nalog na začetku dela je bila izbira časovnega obdobja preučevanja pozidanih površin in temu primerna izbira satelita. Ker smo želeli prikazati prostorsko širjenje pojava, smo izbrali tri časovne prereze: leto 2000, 2013 in 2020. Obdobje 20 let se je zdelo dovolj dolgo, da bodo opazne vidne spremembe, a je bilo treba, zaradi tako širokega razpona, uporabiti tri različne satelite. To so Landsat 5 za leto 2000, Landsat 8 za leto 2013 in Sentinel 2 za leto 2020. V nadaljevanju bomo predstavili nekaj ključnih lastnosti vsakega.
Preglednica 1: Podatki o izbranih satelitskih slikah.
Leto Satelit ID satelitske slike Datum Stopnja
oblačnosti (%) 2000 Landsat 5 LT05_L2SP_190028_20000817_202009
06_02_T1
17. 8. 2000 1 %
2013 Landsat 8 LC81910282013224LGN00 12. 8. 2013 10, 5 %
2020 Sentinel 2 L1C_T33TVM_A019167_20201106T10 0216
6. 11. 2020 2, 2 % Vir: Sentinel Hub, 2021; USGS, 2021.
3. 1. 2. 1. Landsat 5
Satelit Landsat 5 je bil izstreljen leta 1984 in je posredoval podatke o Zemljinem površju slabih 29 let. Okoli Zemlje je krožil na višini 705 km in jo obkrožil vsakih 99 minut, isto lokacijo pa je dosegel na 16 dni. V celotnem obdobju delovanja je Zemljo obkrožil več kot 150.000 krat in svetu posredoval več kot 2,5 milijona posnetkov kopnega. Radiometrična ločljivost satelita je 8 bitna, a so slike v nadaljevanju obdelane in shranjene v 16 bitni ločljivosti. Imel je dva senzorja, in sicer MMS (Multispectral Scanner) in TM (Thematic Mapper). Za analizo smo uporabili senzor TM, ki ima sedem spektralnih pasov z ločljivostjo 30 m, z izjemo termičnega pasu, ki ima ločljivost 120 m, a za nas ni bil pomemben, saj ni bil vključen v analizo. Lastnosti posameznega pasu prikazuje preglednica 2 (Landsat 5, 2021).
9
Preglednica 2: Značilnosti satelita Landsat 5.
Vir: Landsat 5 (TM), 2021.
Podatke satelita smo pridobili s spletne strani United States Geological Survey (USGS), kjer smo se morali predhodno prijaviti. Izbrani posnetki so bili posneti 17. 8. 2000 z 1 % oblačnostjo.
Ker so bili dobljeni v .tiff obliki, jim ni bilo treba spreminjati formata.
3. 1. 2. 2. Landsat 8
Satelit Landsat 8 je bil izstreljen februarja leta 2013 in še vedno deluje. Okoli Zemlje kroži na isti višini kot je satelit Landsat 5, prav tako Zemljo obkroži v 99 minutah in se nad isto lokacijo pojavi na 16 dni. Kljub temu pa na dan zajame več slik, kar povečuje verjetnost zajema brez oblakov. Radiometrična ločljivost je 12 bitna, a so v nadaljevanju slike obdelane, kjer se ločljivost spremeni na 16 bitov (Grant Hamilton, 2018). Ima dva senzorja, to sta OLI (Operational Land Imager) z devetimi spektralnimi pasovi in TIRS (Thermal Infrared Sensor) z dvema spektralnima pasovoma. V analizi smo torej uporabili spektralne pasove senzorja OLI, ki predstavlja izboljšavo preteklim Landsatovim senzorjem, zato lahko pričakujemo nekoliko bolj kakovostne slike, kot pri satelitu Landsat 5. Lastnosti posameznega pasu prikazuje preglednica 3 (Landsat 8 Overview, 2021).
Preglednica 3: Značilnosti satelita Landsat 8.
Kanal Spektralni pas Valovna
dolžina (µm)
Ločljivost (m)
B1 Modra 0.45 – 0.52 30
B2 Zelena 0.52 – 0.60 30
B3 Rdeča 0.63 – 0.69 30
B4 Bližnja infrardeča – NIR 0.76 – 0.90 30 B5 Bližnja infrardeča –
SWIR1
1.55 – 1.75 30
B6 Termična IR 10.40 – 12.50 120
B7 Bližnja infrardeča –
SWIR2 2.08 – 2.35 30
Kanal Spektralni pas Valovna
dolžina (µm)
Ločljivost (m)
B1 Obalni aerosol 0.43 – 0.45 30
B2 Modra 0.45 – 0.51 30
B3 Zelena 0.53 – 0.59 30
B4 Rdeča 0.64 – 0.67 30
B5 Bližnja infrardeča – NIR
0.85 – 0.88 30 B6 Bližnja infrardeča –
SWIR 1
1.57 – 1.65 30 B7 Bližnja infrardeča –
SWIR2
2.11 – 2.29 30
B8 Pankromatsko 0.50 – 0.68 15
B9 Cirus 1.36 – 1.38 30
10 Vir: Landsat 8, 2021.
Posnetki satelita Landsat 8 so bili pridobljeni s spletne strani EO Browser, kjer smo se tudi morali predhodno prijaviti. Najbolj primeren posnetek za izbrano leto 2013 je bil 12. 8. 2013 z 10, 5 % oblačnostjo.
3. 1. 2. 3. Sentinel 2
Satelit Sentinel 2 je bil izstreljen leta 2015 v okviru programa Evropske komisije Copernicus in deluje še danes. Ker gre za najnovejši in najbolj natančen izbrani satelit za našo analizo, smo predpostavili, da bo podal najboljše rezultate. Giblje se na višini 785 m in zajema 13 spektralnih pasov z ločljivostjo od 10 m (vidni spekter in NIR) do 60 m. Radiometrična ločljivost je 12 bitna, a je tako kot pri Landsat 8 spremenjena v 16 bitno (Radiometric resolution, 2021).
Podrobnejše podatke za posamezni spektralni pas prikazuje spodnja preglednica (Sentinel-2, 2021).
Preglednica 4: Značilnosti satelita Sentinel 2.
Vir: Goldblatt, Deininger, Hanson, 2018.
V analizi smo pridobili satelitske slike satelita Sentinel 2 L2A na spletni strani United States Geological Survey (USGS). Najbolj primeren posnetek za leto 2020 je bil 6. 11. 2020, saj je bila oblačnost le 2, 2 %, hkrati pa je posnetek v celoti zajemal naše izbrano območje.
B10 Termični infrardeč senzor 1 – TIRS 1
10.60 – 11.19 100 B11 Termični infrardeč
senzor 2 – TIRS 2
11.50 – 12.51 100
Kanal Spektralni pas Valovna dolžina
(nm) Ločljivost
B1 Aerosoli 443 60
B2 Modra 490 10
B3 Zelena 560 10
B4 Rdeča 665 10
B5 Pas roba rdeče barve
1 – RE 1 705 20
B6 Pas roba rdeče barve
2 – RE 2 740 20
B7 Pas roba rdeče barve 3 – RE 3
783 20
B8 Bližnja infrardeča – NIR
842 10
B8a Pas roba rdeče barve 4 – RE 4
865 20
B9 Vodna para 940 60
B11 Bližnja infrardeča –
SWIR 1 1610 20
B12 Bližnja infrardeča – SWIR 2
2190 20
11
Preden smo satelitske slike za leto 2020 uvozili v program ArcMap, smo jim v programu TerrSet z ukazom GDAL Conversation Utility spremenili zapis iz .jp2 v .rst. To je bil edini ukaz uporabljen v programu TerrSet, zato se bodo vsi ukazi, ki bodo izpostavljeni v nadaljevanju, nanašali na program ArcMap.
Ker ta satelit zajema podatke v določenih pasovih z različno ločljivostjo, je bilo potrebno pred primerjavo z ostalimi slikami satelitov Landsat, spremeniti na skupno prostorsko ločljivost 30 m. To smo storili z orodjem Resample, ki omogoča določanje velikosti rastrskih celic. V nadaljevanju smo tako lahko primerjali zaznane pozidane površine na podlagi posnetkov vseh treh satelitov.
3. 2. Izbrani spektralni indeksi
Med pregledom literature smo našli kar nekaj različnih spektralnih indeksov, ki zaznavajo pozidane površine, a smo na koncu izbrali 4, ki so ustrezali uporabi na razpoložljivih podatkih.
Paziti smo morali, da so bili vsi spektralni pasovi, ki so potrebni za izračun spektralnih indeksov tudi v naših treh izbranih satelitih. Na primer pas TIRS, ki je prisoten pri indeksu EBBI in NDBaI ima v našem primeru samo Landsat 8, zato izračuna ne bi morali narediti na posnetkih Landsata 5 in Sentinela 2. Izbrali smo torej 4 indekse, ki so ustrezali našim izbranim podatkom.
To so Normalized Difference Vegetation Index (v nadaljevanju NDVI), Normalized Difference Built-up Index (v nadaljevanju NDBI), Urban Index (v nadaljevanju UI) in Built-up Index (v nadaljevanju BU). Gre za tri indekse, ki so neposredno namenjeni zaznavanju pozidanih površin (to so NDBI, BU in UI) in eden vegetacijski indeks (to je NDVI). V nadaljevanju bomo vsakega posebej kratko predstavili.
Vse izračune, ki so predstavljeni v nadaljevanju, smo izvedli s pomočjo orodja Raster calculator. Pri vstavljanju formul v orodje je bilo treba paziti, da smo pred vstavljenim spektralnim pasom dodali funkcijo Float, da so bili dobljeni rezultati v decimalnih vrednostih.
Slika 3: Primer vstavljene formule za indeks NDVI v orodju Raster calculator.
3. 2. 1. NDVI indeks
Normalized Difference Vegetation Index (NDVI) je najpogosteje uporabljen vegetacijski indeks. Za njegov izračun se uporabita bližnji infrardeč spektralni pas in pas rdeče vidne svetlobe. Natančneje gre za razmerje med omenjenima spektralnima pasovoma, pri čemer se infrardeča svetloba na zdravi vegetaciji odbija, na ostalih površinah pa ne, zato s tem indeksom lahko zaznamo tudi pozidane površine. Vrednosti indeksa se gibajo med -1 in +1, kjer višje vrednosti kažejo območja z bolj intenzivno vegetacijo, nižje pa pozidano (Kshetri, 2018; Krevs, 2004).
12
Za izračun indeksa NDVI smo torej pri vseh treh satelitih potrebovali pas rdeče vidne svetlobe in bližnji infrardeč pas. V primeru Sentinela 2 sta to pas 4 in pas 8, pri Landsatu 8 pas 4 in pas 5 in pri Landsatu 5 pas 3 in pas 4.
Preglednica 5: Formula NDVI indeksa.
Indeks Formula Sentinel 2 Landsat 8 Landsat 5
NDVI (𝑁𝐼𝑅 − 𝑅𝑒𝑑) (𝑁𝐼𝑅 + 𝑅𝑒𝑑)
(𝑝𝑎𝑠 8 − 𝑝𝑎𝑠 4) (𝑝𝑎𝑠 8 + 𝑝𝑎𝑠 4)
(𝑝𝑎𝑠 5 − 𝑝𝑎𝑠 4) (𝑝𝑎𝑠 5 + 𝑝𝑎𝑠 4)
(𝑝𝑎𝑠 4 − 𝑝𝑎𝑠 3) (𝑝𝑎𝑠 4 + 𝑝𝑎𝑠 3) Vir: Kshetri, 2018.
3. 2. 2. NDBI indeks
Normalized Difference Built-up Index (NDBI) je eden izmed najbolj znanih spektralnih indeksov za pozidane površine. Njegove vrednosti se gibljejo med -1 in +1, pri čemer pozidane površine predstavljajo višje vrednosti (Kshetri, 2018). Za njegov izračun potrebujemo kratkovalovni infrardeči pas SWIR in bližnji infrardeči pas NIR. Pas SWIR ima na urbanih območjih večjo odbojnost kot NIR, zato z njuno kombinacijo lahko zaznamo pozidane površine (Zha, Gao, Ni, 2003).
Za izračun NDBI indeksa smo pri satelitu Sentinel 2 uporabili pas 11 in pas 8, pri Landsat 8 pas 6 in pas 5 in pri Landsatu 5 pas 5 in pas 4.
Preglednica 6: Formula NDBI indeksa.
Indeks Formula Sentinel 2 Landsat 8 Landsat 5
NDBI (𝑆𝑊𝐼𝑅 − 𝑁𝐼𝑅) (𝑆𝑊𝐼𝑅 + 𝑁𝐼𝑅)
(𝑝𝑎𝑠 11 − 𝑝𝑎𝑠 8) (𝑝𝑎𝑠 11 + 𝑝𝑎𝑠 8)
(𝑝𝑎𝑠 6 − 𝑝𝑎𝑠 5) (𝑝𝑎𝑠 6 + 𝑝𝑎𝑠 5)
(𝑝𝑎𝑠 5 − 𝑝𝑎𝑠 4) (𝑝𝑎𝑠 5 + 𝑝𝑎𝑠 4) Vir: Kshetri, 2018.
3. 2. 3. UI indeks
Urban Index (UI) je sestavljen iz srednjega infrardečega pasu SWIR 2 in bližnjega infrardečega pasu NIR. Zajema vrednosti od -1 do +1, kjer višje vrednosti predstavljajo pozidane površine (Setiyono, Hastuti, Widayani, 2017).
V našem primeru smo za izračun UI indeksa iz Santinelovih posnetkov uporabili pas 12 in pas 8, iz Landsata 8 pas 7 in pas 5 ter iz Landsata 5 pas 7 in pas 4.
Preglednica 7: Formula UI indeksa.
Indeks Formula Sentinel 2 Landsat 8 Landsat 5
UI (𝑆𝑊𝐼𝑅 2 − 𝑁𝐼𝑅) (𝑆𝑊𝐼𝑅 2 + 𝑁𝐼𝑅)
(𝑝𝑎𝑠 12 − 𝑝𝑎𝑠 8) (𝑝𝑎𝑠 12 + 𝑝𝑎𝑠 8)
(𝑝𝑎𝑠 7 − 𝑝𝑎𝑠 5) (𝑝𝑎𝑠 7 + 𝑝𝑎𝑠 5)
(𝑝𝑎𝑠 7 − 𝑝𝑎𝑠 4) (𝑝𝑎𝑠 7 + 𝑝𝑎𝑠 4) Vir: Setiyono, Hastuti, Widayani, 2017.
13 3. 2. 4. BU indeks
Built-up Index (BU) indeks je sestavljen iz dveh predhodnih indeksov, in sicer NDBI in NDVI (Kshetri, 2018). Njegove vrednosti se gibljejo med -2 in +2, pri čemer višje vrednosti predstavljajo pozidane površine.
Ker smo za izračun BU indeksa potrebovali že prej omenjene spektralne pasove pri indeksu NDBI in NDVI, smo njune vrednosti le odšteli in dobili rezultat.
Preglednica 8: Formula BU indeksa.
Indeks Formula Sentinel 2 Landsat 8 Landsat 5
BU 𝑁𝐷𝐵𝐼 − 𝑁𝐷𝑉𝐼 (𝑝𝑎𝑠 11 − 𝑝𝑎𝑠 8) (𝑝𝑎𝑠 11 + 𝑝𝑎𝑠 8)
−(𝑝𝑎𝑠 8 − 𝑝𝑎𝑠 4) (𝑝𝑎𝑠 8 + 𝑝𝑎𝑠 4)
(𝑝𝑎𝑠 6 − 𝑝𝑎𝑠 5) (𝑝𝑎𝑠 6 + 𝑝𝑎𝑠 5)
−(𝑝𝑎𝑠 5 − 𝑝𝑎𝑠 4) (𝑝𝑎𝑠 5 + 𝑝𝑎𝑠 4)
(𝑝𝑎𝑠 5 − 𝑝𝑎𝑠 4) (𝑝𝑎𝑠 5 + 𝑝𝑎𝑠 4)
−(𝑝𝑎𝑠 4 − 𝑝𝑎𝑠 3) (𝑝𝑎𝑠 4 + 𝑝𝑎𝑠 3) Vir podatkov: Kshetri, 2018.
3. 3. Ročna metoda nastavljanja meje
Ko smo izračunali spektralne indekse za vsa leta, smo morali dobljenim vrednostim, ki so se gibale približno od -1 do +1, le pri BU indeksu od -2 do +2, nastaviti mejo med pozidanimi in nepozidanimi površinami. Ob pregledu literature smo ugotovili, da te vrednosti niso univerzalne, saj se lahko vrednosti istih površin razlikujejo glede na svetlost satelitske slike, čas snemanja, kakovost slike …
Tako smo meje za vsako sliko nastavili ročno s pomočjo orodja Swipe Layer, tako, da smo v ospredje postavili dobljeno klasifikacijo, v ozadje pa barvni kompozit satelitske slike istega leta, ki nam je kazal »realno stanje«. Kompozite smo naredili z ukazom Composite Bands z vsemi tremi pasovi, ki zajemajo vidno svetlobo (modra, zelena, rdeča). Pri satelitu Sentinel 2 in Landstat 8 so to pasovi 2, 3 in 4, pri satelitu Landsat 5 pa pasovi 1, 2 in 3. Tako smo z drsenjem med obema slojema ugotavljali, ali se meja pozidanih površin na klasifikaciji ujema z barvnim kompozitom. Če se sliki nista ujemali, smo na sloju klasifikacije v zavihku Symbology spreminjali vrednost meje. Na tem mestu smo ugotovili, da nobena mejna vrednost ni ustrezala natančni ločnici med pozidanim in nepozidanim, zato smo mejo nastavili tako, da se je v čim večji meri ujemala z mejami pozidanih površin, ki smo jih zaznali s pomočjo vida na podlagi kompozita. Odločali smo se torej za meje, ki so se nam zdele še dopustne, glede na celotno sliko in ne le na robove posameznih primerov pozidanih površin.
Pri prepoznavanju dejanskega stanja površja bi si lahko, namesto z barvnimi kompoziti satelitskih slik, pomagali z letalskimi ortofoto posnetki, ki imajo boljšo prostorsko ločljivost, vendar se leta posnetkov niso ujemala z izbranimi leti satelitskih posnetkov.
14
Slika 4: Uporaba Swipe Layer orodja.
3. 4. Otsujeva metoda nastavljanja meje
Drugi način nastavljanja meje spektralnih indeksov pa je bil samodejen z uporabo orodja Binary Thresholding Function. Gre za orodje, ki računa na podlagi Otsujeve metode in vrednosti (v našem primeru vrednosti spektralnega indeksa) razdeli v dva razreda. Metoda je bila zasnovana z namenom razlikovanja med ospredjem in ozadjem slike, ki ustvari dva razreda z minimalno varianco znotraj razreda. Eden razred predstavljajo nizke vrednosti, drugi pa visoke (Otsu, 1979).
3. 5. Nadzorovana klasifikacija
Drugi način prepoznave pozidanih površin na satelitskih posnetkih, poleg uporabe spektralnih indeksov je bila nadzorovana klasifikacija. Za bolj učinkovito zaznavo pozidanih površin smo poskusili prepoznati pet kategorij rabe tal: gozd, travniki, njive, voda in pozidano. Za vsako kategorijo smo naredila 20 učnih vzorcev, torej je bilo vseh učnih vzorcev na eni sliki 100. V nadaljevanju smo za te učne vzorce z orodjem Create signatures izdelali spektralne podpise.
Te smo nato uporabili v orodju Maximum Likelihood Classification, ki uporablja metodo največjega verjetja. Gre za eno izmed mnogih metod, ki je pri nadzorovani klasifikaciji pogosto uporabljena in se je pri delu avtorjev Akgun, Eronat in Türkizkazala za najboljšo (2004, str. 5). Pri razvrščanju pikslov v izbrane razrede upošteva povprečja razredov, variance in korelacije med njimi. Tako v prostoru ustvari elipsoidne ploskve enake verjetnosti (Oštir, 2006). Rezultat te metode je bila želena klasifikacija, ki smo jo na koncu še reklasificirali z orodjem Reclassify. V začetku izbrane kategorije gozd, travnik, njive in vodo smo združili v eno kategorijo – nepozidano.
15
Slika 5: Primer učnih vzorcev za leto 2000.
3. 6. Ocena kakovosti klasifikacije
Po izdelavi klasifikacije je bilo treba oceniti njeno natančnost, da smo lahko verodostojno interpretirali rezultate. To smo storili tako, da smo na naključno izbranih točkah na sliki klasifikacije preverili, kakšno je realno stanje in ali se to stanje ujema z razredom klasifikacije.
Realno stanje lahko preverimo na terenu, na letalskih posnetkih ali na drugih natančnejših kartah. Na koncu se glede na število pravilno zaznanih točk izračuna kappa koeficient, ki je merilo za natančnost klasifikacije, kjer je razmernostna napaka zmanjšana za pričakovano naključno ujemanje klasifikacije z dejanskim stanjem (Grlj, 2020). Vrednosti kappa koeficienta se gibljejo med 0 in 1. Nižja kot je vrednost, večja je možnost, da je dobljena klasifikacija rezultat naključja, zato torej strimo k čim višjim vrednostim. Landis in Koch sta kappa koeficient interpretirala po naslednjih merilih (Vujović, 2021).
Preglednica 9: Interpretacija kappa koeficienta.
Kappa koeficient Interpretacija
≤ 0 slabo
0.01 – 0.20 nepomembno
0.21 – 0.40 pošteno
0.41 – 0.60 zmerno
0.61 – 0.80 znatno
0.81 – 1 skoraj popolno
Vir: Vujović, 2021
16
Pred oceno kakovosti klasifikacije je bilo potrebno reklasificirati slike. Površine, ki niso pozidane so dobile novo vrednost 0, površine, ki so pozidane pa 1. To smo storili z orodjem Reclassify. Ko so bile slike reklasificirane, smo z orodjem Create Accuracy Assessment Points vsaki kategoriji (torej 0 in 1) dodelili 125 naključnih točk, za katere smo v nadaljevanju preverili, ali se ujemamo z realnim stanjem rabe tal. Izbrali smo princip vzorčenja "Equalized stratified random", to pomeni enakovredno med dve kategoriji, pozidano in nepozidano.
Slika 6: Točke, ki smo jih uporabili za oceno kakovosti klasifikacije za leto 2000.
Skupno smo torej izdelali 250 točk za posamezno klasifikacijo. V našem primeru, kot že omenjeno, se letalski posnetki časovno niso ujemali z leti, terenski ogled pa se je zdel nesmiseln, zato smo realno stanje razbrali iz izbranih satelitskih posnetkov. Ker smo izdelali 27 klasifikacij, bi postopek pregleda vseh točk na vseh klasifikacijah trajal več ur, zato smo pregledali le realna stanja za vse tri časovne prereze, kar nam je skupno vzelo nekaj manj kot tri ure, in potem te rezultate prenesli na izdelane klasifikacije. Najprej smo morali klasifikacije spremeniti iz rastrske v vektorsko obliko, da smo dobili poligone pozidanih površin, nato pa smo z orodjem Intersect te poligone združili s točkami realnega stanja, ter dobili nov vektorski sloj. Tega smo nato želeli vstavit v orodje Compute Confusion Matrix, ki nam bi samodejno izračunal kappa indeks, a naših novih slojev orodje ni hotelo sprejeti, zato smo uporabili program Microsoft Exel, kamor smo vnesli vse atributivne tabele in s pomočjo funkcij izračunali kappa koeficiente. V nadaljevanju bomo na kratko predstavili postopek izračuna kappa koeficienta s pomočjo Microsoft Exela.
17
Preglednica 10: Uvoz atributivnih tabel v Microsoft Exel.
Sprva smo vse atributivne tabele uvozili v Microsoft Exel (slika 5). V nadaljevanju smo ustvarili nov stolpec (označen z rdečo puščico) in mu določili vrednost glede na razmerje med realnim stanjem ter stanjem zajema klasifikacije. Te rezultate smo nato sešteli in ugotovili koliko površin je bilo zaznanih pravilno (Class=0; GT=0 ali Class=1; GT=1) in koliko narobe (Class=0; GT=1 ali Class=1; GT=0). Za izračun kappa koeficienta smo potrebovali dve novi vrednosti, to sta p_e in p_o, ki smo ju izračunali s formulama na sliki 6 in 7. Vrednost p_o hkrati predstavlja tudi delež pravilno klasificirani točk. Končni kappa koeficient smo dobili po formuli, ki je prikazana na sliki 8.
Preglednica 11: Izračun vrednosti p_e.
18
Preglednica 12: Izračun vrednosti p_o.
Preglednica 13: Izračun kappa indeksa.
3. 7. Primerjava posnetkov po časovnih obdobjih
Ker smo želeli prikazati širjenje pozidanih površin od leta 2000 do 2020, smo dobljene rezultate po časovnih obdobjih na koncu primerjali med seboj.
Zaradi različnih satelitov in njihove razlike v kakovosti posnetkov smo morali pred primerjavo vsem klasifikacijam poenotiti velikost rastrskih celic. Klasifikacije, ki so bile izdelane na podlagi satelitskih posnetkov Landsat 5 in Landsat 8 so imele prostorsko ločljivost 30 x 30
19
metrov, Sentinel 2 pa ponekod tudi 10 x 10 metrov, zato smo Sentinelovim klasifikacijam s pomočjo orodja Resample spremenili velikost rastrskih celic v 30 x 30 metrov.
Za izračun samih površin smo morali klasifikacije spremeniti v vektorsko obliko. To smo storili z orodjem Raster to Polygon in izbrali možnost Simplify Polygons, ki mejo poligonov zračuna z glajenjem rastrskih celic. Na koncu smo v atributivni tabeli vektorskih slojev dodali nov stolpec in z ukazom Area izračunali površine v km2.
3. 8. Uporabljena orodja
Večino analiz smo naredili v programu Arc Map 10. 8. 1. z orodji:
• Resample: določanje velikosti rastrskih celic
• Project: sprememba koordinatnega sistema
• Composite Bands: sestavljanje barvnega kompozita
• Raster calculator: izračun spektralnih indeksov
• Reclassify: reklasificiranje vrednosti
• Swipe Layer: pomoč pri ročnem nastavljanju meje indeksov
• Binary Thresholding Function: orodje za izračun meje po Otsujevi metodi
• Create signatures: izračun spektralnih podpisov
• Maximum Likelihood Classification
• Create Accuracy Assessment Points: izdelava točk za oceno klasifikacije
• Intersect: združevanje točk za oceno klasifikacije in vektorskega sloja klasifikacije
• Compute Confusion Matrix: izračun Kappa indeksa
• Rraster to Polygon: pretvorba rastrskih podatkov v vektorske
• Calculate geometry – area: Izračun površine pozidanih površin Le eno orodje smo uporabili tudi v programu TerrSet:
• GDAL Conversation Utility: pretvorba satelitskih posnetkov iz zapisa .jp2 v .rst
Kot že omenjeno, smo oceno kakovosti klasifikacije s pomočjo preprostih računskih operacij izračunali v programu Microsoft Exel.
20
4. Rezultati
4. 1. Nadzorovana klasifikacija
Metoda nadzorovane klasifikacije nam je, v primerjavi s spektralnimi indeksi, vzela največ časa, saj je bilo treba določiti 100 učnih vzorcev za vsako leto, torej skupno 300 učnih vzorcev.
Največja prednost te metode je, da smo lahko iz pozidanih površin s pomočjo učnih vzorcev izločili vodo, kar ostale klasifikacije z spektralnimi indeksi niso omogočale.
Za boljšo kakovost klasifikacije smo zajeli pet različnih rab tal (gozd, travnik, njive, voda in pozidano) in nato te reklasificirali v dva razreda – pozidano in nepozidano. Vmesna klasifikacija s petimi kategorijami je predstavljena na spodnji karti. Klasifikacija za leto 2000 je vseh pet kategorij lepo zaznala, izstopajo le vode. Vse vodne površine v Ljubljani so sicer pravilno zaznane kot voda, vendar so tudi nekateri drugi deli površja padli pod to kategorijo.
Ta napaka je sicer problematična, a z vidika cilja zaznavanja pozidanih površin to ne vpliva na kakovost naših rezultatov. Enako se je zgodilo tudi pri klasifikaciji za leto 2013, le da je teh napak malo več. Kljub temu pa še vedno ne vplivajo na kakovost končne klasifikacije. Vse ostale kategorije so večinoma pravilno zaznane. Klasifikacija za leto 2020 zgoraj omenjenih težav ni imela in lahko rečemo, da je v primerjavi s prejšnjima letoma najlepše zaznala posamezne kategorije, kar si razlagamo s tem, da gre za bolj kakovostne satelitske posnetke glede prostorske ločljivosti.
Slika 7: Nadzorovana klasifikacija pred reklasifikacijo.
21
Ob podrobnejši analizi reklasificiranih rezultatov po letih smo opazili, da je klasifikacija za leto 2000 dokaj natančna. Slabši kappa koeficient tako ni rezultat napak večjih površin, ampak napak na posameznih pikslih. Klasifikacija za leto 2013 je najboljše zajeta klasifikacija za to leto med vsemi metodami. To nakazuje tudi visok kappa koeficient. Predvsem je razlika v tem, da nadzorovana klasifikacija za leto 2013 ni zaznala ogromnih površin njiv na severnem delu občine pod pozidano, kot so to zaznale vse ostale klasifikacije za to leto. Tudi rezultat za leto 2020 je dokaj dober v primerjavi z ostalimi v nadaljevanju. Kljub temu pa so večje napake na vzhodnem delu občine, na območju Posavskega hribovja, kjer so zaznane večje pozidane površine na območju gozda. Če pogledamo deleže pozidanih površin, opazimo, da delež narašča postopoma z leti, kar odraža logično širjenje mesta.
Po pregledu rezultatov vseh klasifikacij v diplomski nalogi smo ugotovili, da je nadzorovana klasifikacija dala najboljše rezultate. Hkrati je med vsemi izvedenimi metodami le ta dala rezultate, kjer deleži pozidanih površin postopno naraščajo z leti. Kakovost nadzorovane klasifikacije odražajo tudi visoki kappa koeficienti.
Slika 8: Pozidane površine v Mestni občini Ljubljana izračunane z metodo nadzorovane klasifikacije.
22
Preglednica 14: Delež pozidanih površin in kappa koeficient za nadzorovano klasifikacijo.
Nadzorovana
klasifikacija Pozidano (km2) Pozidano (%) Kappa koeficient
2000 56.4273 20,53 % 0.699057
2013 61.98326365000 22.55 % 0.84000000000
2020 67.71806215000 24.64 % 0.74400000000
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
4. 2. Spektralni indeksi z mejo po Otsujevi metodi
Rezultati zaznavanja pozidanih površin z spektralnimi indeksi, ki smo jim nastavili mejo po Otsujevi metodi, že na prvi pogled nekoliko odstopajo od nadzorovanih klasifikacij. Najslabše rezultate smo dobili pri klasifikacijah po Otsujevi metodi iz leta 2020, kjer so pozidane površine zaznane čez celotno občino, brez logične razlage. Težava te metode je, da orodje samodejno razdeli vrednosti na dva dela z mejo glede na razlikovanje ospredja in ozadja slike. Pozidane površine pri vseh spektralnih indeksih, razen pri indeksu NDVI, zajemajo višje vrednosti (okoli 0 do 1), a hkrati lahko te vrednosti zajemajo tudi druge površine, ki bi jih v naši klasifikaciji radi zaznali kot nepozidane. Tako se torej vse te površine z višjimi vrednostmi avtomatično uvrstijo pod pozidane, čeprav niso. Lep primer tega je leto 2020. Hkrati so tudi vse vodne površine pri vseh indeksih zaznane kot pozidane. Ker program samodejno postavi mejo, konkretnih vrednosti ne navajamo. V nadaljevanju bomo natančneje predstavili rezultate po posameznih spektralnih indeksih.
4. 2. 1. NDVI
Klasifikacija za leto 2000 je v pozidane površine dodatno zajela še nekatere njive in osojna pobočja hribov, enako tudi pri letu 2013, vendar v večji meri. Tu lahko zaznamo, da je celotno Ljubljansko polje zajeto pod pozidano. To si razlagamo s tem, da so bile njive tisti dan zelo osušene, zato so se njihove spektralne vrednosti približevali bolj pozidanim površinam, kot zeleni vegetaciji. Tudi v nadaljevanju pri drugih indeksih smo imeli s tem težave. Ob podrobnejšem pregledu smo ugotovili, da so se na klasifikaciji za leto 2013 pod nepozidano uvrstila le območja, ki so bila zelena, torej travniki in gozdovi. Leto 2020 je med vsemi indeksi z Otsujevo metodo najboljše zajeto z NDVI indeksom, kljub temu da so pod pozidano zajeti tudi nekateri travniki in njive.
Kappa indeksi je pri letu 2013 najvišji, čeprav že na prvi pogled daje najslabši rezultat med tremi leti. To si razlagamo s tem, da kljub temu da je klasifikacija zajela preveč pozidanih površin, so bile vse tiste, ki so pozidane, pravilno klasificirane, torej so napake nastale le pri nepozidanih, a tudi to ne pri vseh.
Ob pregledu vseh indeksov z mejo po Otsujevi metodi, smo ugotovili, da je indeks NDVI v celoti dal najboljše rezultate. Če bi odstranili zgoraj omenjene napake, bi lahko tudi delež pozidanih površin po letih smiselno naraščal.
23
Slika 9: Pozidane površine v Mestni občini Ljubljana izračunane z NDVI indeksom in mejo nastavljeno po Otsujevi metodi.
Preglednica 15: Delež pozidanih površin in kappa koeficient za indeks NDVI in mejo nastavljeno po Otsujevi metodi.
NDVI Pozidano (km2) Pozidano (%) Kappa koeficient
2000 70.81974 25.76841 0.670554
2013 105.3327 38.32462 0.722134
2020 59.12441 21.51638 0.624875
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
4. 2. 2. NDBI
Indeks NDBI je za razliko od indeksa NDVI slabše zajel pozidane površine. Leto 2000 je zajelo še več njivskih površin pod pozidano, 2013 pa je vse površine, tudi travnike zajel pod pozidane.
Izrazito so torej ostali nepozidani le gozdovi. Leto 2020 je najslabše zajeto leto med vsemi metodami. Klasifikacija je nesmiselna, saj je več kot polovica površin v občini zajeta pod pozidano. Temu primerna je tudi kakovost klasifikacije izražena s kappa koeficientom, ki znaša 0,18, kar je najslabše med vsemi klasifikacijami v tej diplomski nalogi.
Delež pozidanih površin sicer čez leta narašča, a opazno narašča prehitro in ne realno.
24
Slika 10: Pozidane površine v Mestni občini Ljubljana izračunane z NDBI indeksom in mejo nastavljeno po Otsujevi metodi.
Preglednica 16: Delež pozidanih površin in kappa koeficient za indeks NDBI in mejo nastavljeno po Otsujevi metodi.
NDBI Pozidano (km2) Pozidano (%) Kappa koeficient
2000 90.6262 32.97544 0.584
2013 123.3897 44.89707 0.668414
2020 155.0159 56.41967 0.184083
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
4. 2. 3. UI
UI indeks je za leto 2000 kar dobro zajel pozidane površine, predvsem v vzhodnem delu občine, na območju Posavskega hribovja. Kljub temu da je zajetih tudi nekaj njivskih površin, bi lahko ta klasifikacija odražala dokaj realno stanje. Vendar pa sta v nadaljevanju ostali dve leti precej slabše zajeti. Napake so podobne kot pri indeksu NDBI, le da je natančnost vseeno malce boljša.
Delež pozidanih površin tu ponovno smiselno narašča, a to ni odraz realnega stanja, ampak večjih napak.
25
Slika 11: Pozidane površine v Mestni občini Ljubljana izračunane z UI indeksom in mejo nastavljeno po Otsujevi metodi.
Preglednica 17: Delež pozidanih površin in kappa koeficient za indeks UI in mejo nastavljeno po Otsujevi metodi.
UI Pozidano (km2) Pozidano (%) Kappa koeficient
2000 76.30826 27.76518 0.676019
2013 116.165 42.26652 0.715207
2020 130.9074 47.64662 0.456282
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
4. 2. 4. BU
BU indeks daje podobne rezultate kot UI in tudi podobno natančnost. Pri letu 2000 je ponovno lepo zajeto območje vzhodnega dela občine v Posavskem hribovju, kar pri ostalih dveh letih ni.
Leto 2013 je dodatno zajelo še vse travnate površine pod pozidano na tem območju, medtem ko je leto 2020 zajelo večji del tega hribovja pod pozidano, brez logične razlage.
Delež pozidanih površin ponovno narašča, vendar zaradi napak.
26
Slika 12: Pozidane površine v Mestni občini Ljubljana izračunane z BU indeksom in mejo nastavljeno po Otsujevi metodi.
Preglednica 18: Delež pozidanih površin in kappa koeficient za indeks BU in mejo nastavljeno po Otsujevi metodi.
BU Pozidano (km2) Pozidano (%) Kappa koeficient
2000 77.02032 28.02326 0.651359
2013 116.7675 42.48548 0.714914
2020 128.3766 46.72568 0.504887
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
4. 3. Spektralni indeksi z ročno nastavljeno mejo
Rezultati spektralnih indeksov z ročno nastavljeno mejo, se veliko bolj smiselni. Težava te metode je, da smo težko določili neko optimalno mejo, za katero bi lahko potrdili, da je edino pravilna. V tem primeru je meja zelo odvisna od posameznika, ki le-to nastavlja. Vrednosti spektralnega indeksa imajo do devet decimalnih vrednosti, kar pomeni, da že majhen premik meje spremeni končni rezultat klasifikacije.
27
Ker imajo satelitske slike različne lastnosti, so meje tudi znotraj istega indeksa v različnih letih postavljene na različnih vrednostih, torej ne gre za neko poenoteno metodo. Mejo smo nastavljali po principu kar se da lepega zajema pozidanih površin, hkrati pa čim manjših napak predvsem v zvezi z njivskimi površinami.
Graf 1: Ročno nastavljene mejne vrednosti med pozidanimi in nepozidanimi površinami na klasifikacijah s spektralnimi indeksi.
Ta metoda nam daje nekakšno dvojnost rezultatov. Po eni strani je natančnost klasifikacij nekoliko slabša kot pri prejšnji metodi, vendar če poznamo satelitske slike Mestne občine Ljubljana, vemo, da so te klasifikacije veliko bolj natančne. Podobno napako smo omenili že pri Otsujevi metodi. Natančnost je na klasifikacijah, ki zajamejo preveč pozidanih površin, lahko boljša, ker so točke, kjer je realno stanje res pozidano, vse pravilno klasificirane, odstopajo le tiste, ki so zaznane kot pozidano, pa realno stanje ni pozidano. V takih primerih bi načeloma morali izbrati večje število točk za oceno klasifikacije, kjer bi več točk lahko padlo na narobe zaznana območja. V nadaljevanju bomo podrobneje predstavili rezultate po spektralnih indeksih.
28 4. 3. 1. NDVI
Indeks NDVI se je podal dokaj smiselne rezultate, izstopa le leto 2020, kjer so bile vrednosti pozidanih površin in njiv zelo blizu. Tako smo zaradi želje po čim manjšem zajemu njiv mejo nastavili tako, da smo odstranili tudi nekaj pozidanih površin v kompaktnem središču mesta.
To odraža tudi delež pozidanih površin za leto 2020, ki je manjši od leta 2000 in 2013. Hkrati je leto 2013 ponovno zajelo del Ljubljanskega polja pod pozidano, kar pa je še vedno bolje, kot pri spektralnih indeksih z mejo po Otsujevi metodi. Vsa tri leta imajo pomanjkljivo zaznane pozidane površine naselji v Posavskem hribovju. Razlog je najverjetneje ta, da so naselja dokaj majhna in obdana z vegetacijo, hkrati pa je prostorska ločljivost posnetkov 30 x 30 m, zato se posamezne hiše hitreje zlijejo z okolico, ki ni pozidana.
Slika 13: Pozidane površine v Mestni občini Ljubljana izračunane z NDVI indeksom in ročno nastavljeno mejo.
Preglednica 19: Delež pozidanih površin in kappa koeficient za indeks NDVI in ročno nastavljeno mejo.
NDVI Pozidano (km2) Pozidano (%) Kappa koeficient
2000 45.40582 16.51957 0.522673
2013 56.92687 20.71076 0.683955
2020 43.32358 15.76631 0.4908
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
29 4. 3. 2. NDBI
Leto 2000 daje pri NDBI indeksu najslabše rezultate. Pozidane površine so zajete pomanjkljivo, medtem ko se v tej kategoriji že pojavljajo večje površine njiv. Slabo stanje odraža tudi kappa koeficient z vrednostjo 0.35. Na drugi strani pa ima leto 2013 eno izmed najvišjih vrednosti kappa koeficienta, kar je zelo vprašljivo. Že na prvi pogled opazimo, da je leto 2013 ponovno zajelo preveč pozidanih površin na Ljubljanskem polju, zato razlog za dobro natančnost pripisujemo ponovno temu, da so vse pozidane površine zajete pod kategorijo pozidano. Leto 2020 spada pod eno izmed najslabše zajetih klasifikacij. Mejo bi sicer lahko prestavili nižje in zajeli več pozidanih površin, a bi s tem pod pozidane pridobili tudi večino ostale občine.
Rezultat za leto 2020 je torej odraz neke optimalne meje, ki zajame pozidne površine pravilno le v središču mesta, ostali deli pa so večinoma zaznani narobe.
Slika 14: Pozidane površine v Mestni občini Ljubljana izračunane z NDBI indeksom in ročno nastavljeno mejo.
Preglednica 20: Delež pozidanih površin in kappa koeficient za indeks NDBI in ročno nastavljeno mejo.
NDBI Pozidano (km2) Pozidano (%) Kappa koeficient
2000 41.652 15.15535 0.357216
2013 78.69304 28.63162 0.565329
2020 22.36632 8.140168 0.222905
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
30 4. 3. 3. UI
UI indeks je dal nekoliko boljše rezultate kot NDBI, predvsem za leto 2020. Leto 2000 ima podobne rezultate kot indeks NDBI, medtem ko je leto 2013 malo boljše zajelo Posavsko hribovje in nekatere njive pod nepozidane površine.
Slika 15: Pozidane površine v Mestni občini Ljubljana izračunane z UI indeksom in ročno nastavljeno mejo.
Preglednica 21: Delež pozidanih površin in kappa koeficient za indeks UI in ročno nastavljeno mejo.
UI Pozidano (km2) Pozidano (%) Kappa koeficient
2000 39.60237 14.40854 0.406033
2013 56.03571 20.38711 0.601678
2020 24.791 9.022221 0.343787
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
31 4. 3. 4. BU
Izmed vseh spektralnih indeksov, ki smo jim mejo nastavljali ročno, smo mnenja, da je BU indeks dal najboljši približek realnemu stanju. Za leto 2000 nam je celo uspelo izločiti manjši del Ljubljanice iz pozidanih površin, kar se pri ostalih spektralniih indeksih ni dalo, saj so bile vrednosti vodnih površin zelo podobne vrednostim pozidanih. Edina težava pri letu 2000 je to, da klasifikacija ni zajela pozidanih površin v Posavskem hribovju. Leto 2013 je ponovno kot vse klasifikacije, zajel njive na Ljubljanskem polju pod pozidano. Ker se je to zgodilo večinoma na vseh klasifikacijah, je razlog najverjetneje v tem, da so bile njive tisti dan močno izsušene.
To lahko sklepamo tudi po tem, da na Ljubljanskem barju z njivami nismo imeli toliko težav, ker to območje velja za bolj mokrotno. S tem se vegetacija ne izsuši in ohranja svojo pravo barvo, ki jo zazna satelit. Klasifikacija za leto 2020 je zajeta mnogo bolje kot pri ostalih spektralnih indeksih. Sicer je na Ljubljanskem barju pod pozidano zajetih tudi nekaj več njiv, a glede na prejšnje rezultate to ni tako veliko odstopanje. Če bi ostranili napake pri letu 2013, bi lahko dobili dokaj smiselno časovno zaporedje večanja pozidanih površin v Mestni občini Ljubljana.
Slika 16: Pozidane površine v Mestni občini Ljubljana izračunane z BU indeksom in ročno nastavljeno mejo.
32
Preglednica 22: Delež pozidanih površin in kappa koeficient za indeks BU in ročno nastavljeno mejo.
BU Pozidano (km2) Pozidano (%) Kappa koeficient
2000 37.76438 13.73981 0.406033
2013 62.63267 22.78703 0.618333
2020 38.40884 13.9784 0.420705
Vir podatkov: Sentinel Hub, 2021; USGS, 2021.
33
5. Sklep
V diplomski nalogi smo želeli s pomočjo preprostih klasifikacijskih postopkov prikazati širjenje pozidanih površin med leti 2000, 2013 in 2020 v Mestni občini Ljubljana in primerjati uspešnost različnih metod. Prvi cilj naloge je bil poiskati primerne satelitske slike za izbrana leta. Ker smo v analizo zajeli tri časovne prereze, smo morali uporabiti satelitske slike treh različnih satelitov. Za leto 2000 smo izbrali slike satelita Landsat 5, za leto 2013 slike satelita Landsat 8 in za leto 2020 slike satelita Sentinel 2. Že zaradi različnih starosti satelitskih sistemov smo lahko opazili razlike. Sentinel 2 je zaradi boljše prostorske ločljivosti dal bolj natančne rezultate z vidika kakovosti zaznavanja pozidanih površin, kar smo opazili pri zajemu pozidanih površin znotraj mestnega središča. Ta je namreč zaznal tudi mestne parke kot zelene površine in jih zajel v nepozidano, medtem ko so klasifikacije za leto 2000 in 2013 mestno središče zajele kot homogeno enoto pozidanih površin. Drugi cilj se je navezoval na izbiro primerne nadzorovane klasifikacije in spektralnih indeksov. Za nadzorovano klasifikacijo smo izbrali metodo največje verjetja, ki se je v predhodnih študijah izkazala kot dobra, med spektralnimi indeksi pa smo izbrali indeks NDVI, NDBI, UI in BU. Prva metoda zajema pozidanih površin je bila nadzorovana klasifikacija, pri kateri smo kot prednost prepoznali enostaven zajem vodnih površin pod kategorijo nepozidano. Končni rezultati te metode so bili izmed vseh metod najboljši in hkrati kažejo logično naraščanje pozidanih površin skozi leta.
Druga metoda zajema pozidanih površin je bila uporaba štirih spektralnih indeksov, ki smo jim meje med pozidanimi in nepozidanimi površinami nastavili avtomatično s pomočjo Otsujeve metode. Metoda vse vrednosti spektralnega indeksa razdeli v dva razreda glede na vrednost piksla. Višje vrednosti so predstavljale pozidane površine, nižje pa nepozidane, z izjemo indeksa NDVI, ki deluje ravno obratno. Metoda se ni izkazala za najboljšo, saj smo ugotovili, da visokih vrednosti nimajo le površine, ki so pozidane, ampak tudi druge. Predvsem voda in njive. Če pa bi morali izpostaviti spektralni indeks z najboljšimi rezultati te metode, bi izpostavili NDVI. Pri tretji metodi smo istim indeksom mejo nastavili ročno s pomočjo podrobnejšega pregleda izdelane klasifikacije. S tem smo lahko iz skupine pozidano dodatno izločili večje njivske površine, ki so uvrstile v kategorijo pozidano, čeprav to v celoti ni uspelo.
Pomanjkljivost te metode je, da vodnih površin nismo morali uvrstiti pod nepozidano. Najbolj smiselne rezultate te metode je dal indeks BU, kar je pričakovano, glede na to, da je med vsemi štirimi najbolj kompleksen. Presenetil je tudi indeks NDVI, ki je v osnovi indeks za zaznavanje zdrave vegetacije. Najverjetneje je to tudi razlog, saj gre za velik kontrast med vegetacijo in pozidanimi površinami. Z Otsujevo metodo in ročno nastavljenimi mejami smo torej izpolnili tretji cilj.
Pri vseh metodah pa smo se srečali z nekaterimi podobnimi težavami in napakami. Zaradi srednje prostorske ločljivosti satelitskih posnetkov je klasifikacija težje zajela samotne hiše ali manjša naselja. To je opazno na območju Posavskega hribovja na vzhodu delu občine, kjer so pozidane površine slabo zaznane. Večje težave so povzročale tudi njive v okolici Ljubljane, ki so se uvrščale med kategorijo pozidano. Razlog za to je lahko ta, da njive v tistem času niso bile obdelane in se s tem bolj približale značilnostim surove, pozidane površine, ali pa so bile zaradi vremenskih razmer izsušene. Največ težav z njivami smo imeli pri satelitskih posnetkih iz leta 2013. Zanimivo dejstvo je tudi to, da so vrednosti kakovosti klasifikacij, izražene s kappa keoficientom, pri nekaterih klasifikacijah dokaj visoke, čeprav lahko že na prvi pogled opazimo, da je zajetih preveč pozidanih površin. Razlog za višjo vrednost kappa keoficienta je verjetno v tem, da so točke za ocenjevanje kakovosti klasifikacije po celotnem preučevanem območju porazdeljene enakomerno, medtem ko so napake zgoščene na enem delu, v največ primerih na območju Ljubljanskega polja. Poleg tega so pri klasifikacijah, ki so zajele preveč pozidanih površin, vse točke, na katerih je tudi realno stanje pozidano, kvalificirano pravilno.
34
Zaradi teh napak pri oceni se torej nismo mogli sklicevati le na ta rezultat. Četrti cilj, da preverimo kakovost klasifikacij s kappa koeficientom je torej bil izveden, a z rezultati nismo bili najbolj zadovoljni.
Zadnji cilj je bil ključni cilj celotne naloge in sicer v obliki kart prikazati dobljene rezultate in jih interpretirati. Prvi del cilja nam ni delal težav, medtem ko je bila interpretacija nekoliko bolj zahtevna, saj je zahtevala natančen vizualni pregled dobljenih klasifikacij in nato pomislek čemu tak rezultat.
V prihodnjih analizah daljinskega zaznavanja pozidanih površin bi se raje odločili za preučevanje satelitskih slik enega ali največ dveh satelitov, saj se le-te med seboj razlikujejo in dajejo različno kakovostne slike, kar se odraža tudi pri končnih rezultatih. Če že uporabljamo slike različnih satelitov, bi v prihodnje pazili na to, da imajo sateliti čim več istih spektralnih pasov snemanja, kar nam omogoča večjo svobodo pri izbiri spektralnih indeksov za pozidane površine, saj za določene indekse potrebujemo spektralne pasove, ki jih nekateri sateliti nimajo.
Poleg tega se je smiselno v prihodnjih analizah bolje prepričati, v kakšnem vremenskem stanju so izbrane satelitske slike posnete. Ob daljših sušah ali po drugi strani po daljših deževnih obdobjih imajo lahko iste površine različne spektralne odboje, kar v nadaljevanju pomeni različne klasifikacije. Glede težav z vodnimi površinami, ki so se skoraj pri vseh klasifikacijah uvrstile med pozidano, pa bi v prihodnje lahko z spektralnim indeksom za vodne površine le-te izločili.
Pred začetkom dela smo si zadali štiri hipoteze, ki smo jih nato z dobljenimi rezultati potrdili ali zavrgli. Nekatere klasifikacije, izračunane s spektralnimi indeksi, niso konkretno nakazovale, da so se pozidane površine v preučevanem obdobju v Mestni občini Ljubljana povečale. Vzrok za to so napake v zaznavi površja. Tako so imele v nekaterih primerih klasifikacije za leto 2000 več pozidanih površin kot za leto 2020. Prvo hipotezo smo torej zavrnili. Drugo hipotezo, da bo metoda zaznavanja pozidanih površin s spektralnimi indeksi dala boljše rezultate kot nadzorovana klasifikacija, smo zavrgli. Metoda nadzorovane klasifikacije je namreč edina metoda, ki je v celoti zajela vodne površine pod nepozidane, medtem ko s spektralnimi indeksi to ni bilo mogoče, oz. nam je to uspelo le pri indeksu BU z ročno nastavljeno mejo. Hkrati je tudi kappa koeficient pri nadzorovani klasifikaciji za vsa tri leta, glede na ostale klasifikacije, boljši. Tretjo hipotezo, da bo indeksu BU zaradi svoje kompleksnosti, podal boljše rezultate, smo potrdili. Zadnjo hipotezo, da bo ocena kakovosti klasifikacije po kappa koeficientu vsaj pri eni klasifikaciji nad 0.8 smo potrdili, saj je kakovost pri nadzorovani klasifikaciji za leto 2013 presegla zadano vrednost.
35
5. Summary
In this paper we wanted to use simple classification procedures to show the expansion of built- up areas between years 2000, 2013 and 2020 in the Municipality of Ljubljana and to compare different methods. The first goal was to find suitable satellite images for selected years. Since we included three time sections in the analysis, we had to use satellite images of three different satellites. For 2000 we selected images of the Landsat 5 satellite, for 2013 images of the Landsat 8 satellite and for 2020 images of the Sentinel 2 satellite. The second goal was related to the selection of an appropriate supervised classification and spectral indices. For the supervised classification, we chose the maximum likelihood method, which proved to be good in previous studies, and among the spectral indices, we chose the NDVI, NDBI, UI, and BU index. The first method of capturing built-up areas was the supervised classification, in which we recognized the simple capture of water surfaces under the category of unbuilt-up as an advantage. The final results of this method were the best of all methods and at the same time show a logical increase in built-up areas over the years. Another method of capturing built-up areas was the use of four spectral indices, for which the boundaries between built-up and not built-up areas were set automatically using the Otsu method. The method did not prove to be the best, as we found that not only the surfaces that are built in, but also others have high values.
Especially water and fields. However, if we were to expose the spectral index with the best results of this method, we would expose the NDVI. In the third method, we set the limit for the same indexes manually using a more detailed review of the produced classification. With this, we were able to additionally exclude larger arable land from the group of built-up areas, which were classified in the category of built-up areas, although this was not entirely successful. The disadvantage of this method is that we did not have to classify water surfaces as unbuilt. The most meaningful results of this method were given by the BU index, which expected, given that it is the most complex of them all. We were also surprised by the NDVI index, which is basically an index for the perception of healthy vegetation. This is most likely also the reason, as there is a great contrast between the vegetation and the built-up areas. With the Otsu method and manually set limits, we thus fulfilled the third goal. A fourth goal to check the quality of the kappa coefficient classifications was implemented, but we were not the most satisfied with the results. The quality values of the classifications, expressed by the kappa coefficient, are quite high in some classifications, although at first glance we can see that too many built-up areas are covered. The last goal was the main goal of the whole task, namely to show the obtained results in the form of maps and to interpret them. The first part of the goal did not cause us any problems, while the interpretation was a bit more demanding, as it required a detailed visual review of the obtained classifications and then the thought of why such a result.
Before starting the work, we set ourselves four hypotheses, which we then confirmed or rejected with the obtained results. Some classifications, calculated with spectral indices, did not specifically indicate that the built-up areas in the City of Ljubljana increased during the study period. This is caused by errors in surface perception. Thus, in some cases, the classifications for 2000 had more built-up areas than for 2020. The first hypothesis was therefore rejected. The second hypothesis that the method of detecting built - up areas with spectral indices will give better results than the supervised classification was rejected. Namely, the method of supervised classification is the only method that completely covered water surfaces under unbuilt, while with spectral indices this was not possible, or we only managed this with the BU index with a manually set limit. At the same time, the kappa coefficient in the supervised classification for all three years, compared to other classifications, is better. We confirmed the third hypothesis that the BU index will give better results due to its complexity. The last hypothesis that the assessment of the quality of the classification according to the kappa coefficient in at least one
