Reflektancia vizsgálata

Töltsük be a következő rétegeket:

  • Sentinel-2 műholdfelvételek Tihanyra (RGB, MS), 2019
  • Corine vektoros felszínborítás Tihanyra, 2018

1. ábra: Sentinel-2 RGB színes műholdfelvétel, Tihany, 2019.
1. ábra: Sentinel-2 RGB színes műholdfelvétel, Tihany, 2019.

2. ábra: Corine Land Cover felszínborítás, Tihany, 2018.
2. ábra: Corine Land Cover felszínborítás, Tihany, 2018.

A vetületek harmonizálását követően - jelen esetben a vektor réteget transzformáltuk a raszter réteg vetületéhez -, tisztítsuk meg a réteget az esetleges forgácspoligonoktól, amelyek a településhatárhoz igazított földhasználatból következhetnek.

Ezt követően vonjuk össze multipoligonná az azonos kóddal rendelkező elemeket (Geoprocesszing eszközök / Összevon).

Írjunk egy kódot a rétegek betöltésére, illetve a zóna statisztikák számítására.

path = "d:/indikatrix-org/reflektancia/" 

raster_filename = "sentinel2_tihany_rgb"
vector_filename = "clc_2018_tihany_utm_multi"

rasterlayer = QgsRasterLayer(path + raster_filename + ".tif", raster_filename)
vectorlayer = QgsVectorLayer(path + vector_filename + ".shp", vector_filename, "ogr")

#processing.algorithmHelp("qgis:zonalstatistics")

parameters = {'INPUT_RASTER': rasterlayer,
              'RASTER_BAND': 1,
              'INPUT_VECTOR': vectorlayer,
              'COLUMN_PREFIX': '1_',
              'STATS': "2, 4"
              }

processing.run("qgis:zonalstatistics", parameters)

QgsProject.instance().addMapLayer(vectorlayer)
del(vectorlayer)

Bővítsük úgy a kódot, hogy az valamennyi raszter sávon készítsen statisztikákat (2: átlag, 4: szórás):

for band in range (1, rasterlayer.bandCount() + 1): 
    #processing.algorithmHelp("qgis:zonalstatistics")
    parameters = {'INPUT_RASTER': rasterlayer,
                  'RASTER_BAND': band,
                  'INPUT_VECTOR': vectorlayer,
                  'COLUMN_PREFIX': str(band) + "_",
                  'STATS': "2, 4"
                  }
    processing.run("qgis:zonalstatistics", parameters)

Bővítsuk tovább a kódot úgy, hogy végigfusson valamennyi felszínborítás osztályon, illetve valamennyi sávon (RGB esetén 3...). Készítsünk ábrát az eredményeinkről, illetve fordítsuk magyarra a jelmagyarázatot:

from qgis.core import *
import processing

path = "d:/indikatrix-org/reflektancia/" 

raster_filename = "sentinel2_tihany_rgb"

vector_filename = "clc_2018_tihany_utm_multi"

rasterlayer = QgsRasterLayer(path + raster_filename + ".tif", raster_filename)
vectorlayer = QgsVectorLayer(path + vector_filename + ".shp", vector_filename, "ogr")

vectorlayer.loadNamedStyle("d:/indikatrix-org/reflektancia/clc_jelrendszer_vektor_hu.qml")
vectorlayer.triggerRepaint()
renderer = vectorlayer.renderer()
if renderer.type() == "categorizedSymbol":
    corine_dict = {}
    for cat in renderer.categories():
        cat_split = cat.label().split(" - ")
        corine_dict[cat_split[0]] = cat_split[1]

for band in range (1, rasterlayer.bandCount() + 1): 
    #processing.algorithmHelp("qgis:zonalstatistics")
    parameters = {'INPUT_RASTER': rasterlayer,
                  'RASTER_BAND': band,
                  'INPUT_VECTOR': vectorlayer,
                  'COLUMN_PREFIX': str(band) + "_",
                  'STATS': "2, 4"
                  }
    processing.run("qgis:zonalstatistics", parameters)


import matplotlib.pyplot as plt

fig, ax = plt.subplots()

### here comes the nanometers of bands (BGR)
x = [665, 560, 490]

select_field = "clc2018"
idx_label = vectorlayer.fields().lookupField(select_field)

corine_classes = [feature for feature in vectorlayer.getFeatures()] 

for corine_class in corine_classes:
    y = []
    color = next(ax._get_lines.prop_cycler)['color']
    for band in range (1, rasterlayer.bandCount() + 1): 
        statid = str(band) + "_mean"
        idx = vectorlayer.fields().lookupField(statid)
        y.append(corine_class.attributes()[idx])
    #classlabel = corine_class.attributes()[idx_label]
    classlabel = corine_dict[str(corine_class.attributes()[idx_label])]
    line, = ax.plot(x, y, label = str(classlabel), color = color)
    # in case of adding limits by stdev
    #line.set_dashes([2, 2, 10, 2])  # 2pt line, 2pt break, 10pt line, 2pt break
    
ax.legend(loc = 'upper right', fontsize = 'xx-small')
plt.show()

Az osztálynevek magyarra fordítását egy a rétegre ráolvasott stílusfájl segítségével végeztük, a 'corine_dict' asszociatív tömb változóba (angolul map vagy dictionary) gyűjtöttük össze a kódokhoz tartozó feliratokat.

A készített ábráról jól leolvasható, hogy például a vízfelszíneknek a legnagyobb a visszaverődési értéke a kék színtartományban:

3. ábra: Felszínborítási osztályok visszaverődésének vizsgálata színes RGB Sentinel-2 műholdfelvételen.
3. ábra: Felszínborítási osztályok visszaverődésének vizsgálata Sentinel-2 műholdfelvételen.

A kód finomhangolásával könnyen elkészíthetjük a visszaverődési ábra multispektrális változtatát:

raster_filename = "sentinel2_tihany_msi"
x = [443, 490, 560, 665, 705, 740, 783, 842, 865, 945, 1375, 1610, 2190]

4. ábra: Felszínborítási osztályok visszaverődésének vizsgálata multispektális (13 sáv) Sentinel-2 műholdfelvételen.
4. ábra: Felszínborítási osztályok visszaverődésének vizsgálata multispektális (13 sáv) Sentinel-2 műholdfelvételen.

Módosítsuk úgy a kódot, hogy az ábrán csak bizonyos osztályok szerepeljenek:

class_selector = (512, 231, 311)
select_expr = select_field + " IN " + str(class_selector)

request = QgsFeatureRequest().setFilterExpression(select_expr)

corine_classes = [feature for feature in vectorlayer.getFeatures(request)]
#corine_classes = [feature for feature in vectorlayer.getFeatures()]

5. ábra: Kiválasztott felszínborítási osztályok - tó, erdő, legelő - visszaverődésének vizsgálata színes Sentinel-2 műholdfelvételen.
5. ábra: Kiválasztott felszínborítási osztályok - tó, erdő, legelő - visszaverődésének vizsgálata színes Sentinel-2 műholdfelvételen.

Itt még nem ér véget a történet, - bár a diszciplína számos cikkéből elhagyják a szórást, azért érdemes megvizsgálni. - Akár kevesebb felszínborítási osztályon.

Normális eloszlást feltételezve 1 szigmás határon belül szerepel az átlag megállapításához felhasznált adatok 68 %-a, 2 szigma esetén ez már kicsit több mint 95 %.

6. ábra: Kiválasztott felszínborítási osztályok - tó, erdő, legelő - visszaverődésének vizsgálata színes Sentinel-2 műholdfelvételen, 68 %-os szórásintervallumokkal.
6. ábra: Kiválasztott felszínborítási osztályok - tó, erdő, legelő - visszaverődésének vizsgálata színes Sentinel-2 műholdfelvételen, 68 %-os szórásintervallumokkal.

A fenti ábra elkészítéséhez módosítsuk a kódunkat a következőképp:

    for band in range (1, rasterlayer.bandCount() + 1): 
        statid = str(band) + "_mean"
        idx_mean = vectorlayer.fields().lookupField(statid)  
        y.append(corine_class.attributes()[idx_mean])
        
        statid = str(band) + "_stdev"
        idx_stdev = vectorlayer.fields().lookupField(statid)        
        lci.append(corine_class.attributes()[idx_mean] - 1 * corine_class.attributes()[idx_stdev])
        uci.append(corine_class.attributes()[idx_mean] + 1 * corine_class.attributes()[idx_stdev])
        
    #classlabel = corine_class.attributes()[idx_label]
    classlabel = corine_dict[str(corine_class.attributes()[idx_label])]
    line, = ax.plot(x, y, label = str(classlabel), color = color)
    # in case of adding limits by stdev, plot them by dash line
    uci_line, = ax.plot(x, uci, color = color)
    uci_line.set_dashes([2, 2, 10, 2])  # 2pt line, 2pt break, 10pt line, 2pt break
    lci_line, = ax.plot(x, lci, color = color)
    lci_line.set_dashes([2, 2, 10, 2])  

7. ábra: Kiválasztott felszínborítási osztályok - tó, erdő, legelő - visszaverődésének vizsgálata színes Sentinel-2 műholdfelvételen, 95 %-os szórásintervallumokkal.
7. ábra: Kiválasztott felszínborítási osztályok - tó, erdő, legelő - visszaverődésének vizsgálata színes Sentinel-2 műholdfelvételen, 95 %-os szórásintervallumokkal.

A 7. ábrán jól látható, hogy a határok meglehetősen összemosódtak 95 %-os szórásintervallumok esetén.

© Dr. Wirth Ervin