153ZODH / 15. cvičení
Objektově orientovaná klasifikace
Osnova
Toto cvičení je zaměřeno na možnosti objektově orientované klasifikace družicových dat. V tomto ohledu navazuje na cvičení věnované klasifikaci, řizené klasifikaci a především segmentaci obrazu.
Seznam příkazů
- i.group
- g.region
- i.segment
- r.to.vect
- v.db.dropcolumn
- v.db.addcolumn
- v.db.update
- v.rast.stats
- r.texture
- v.to.db
- v.db.select
- v.extract
- v.class.mlpy
- v.to.rast
Příprava
Nejprve nastavíme výpočetní region (g.region). Vzhledem k tomu, že je proces objektově orientované klasifikace poměrně výpočetně náročný, zvolíme menší region. V našem případě předpokládáme existenci již definované regionu uloženého jako 'ceske-budejovice-10km' (jde o region s offsetem 10km kolem města České Budějovice, viz návod).
g.region region=ceske-budejovice-10km -p
projection: 1 (UTM) zone: 33 datum: wgs84 ellipsoid: wgs84 north: 5439839.30457264 south: 5410176.48726338 west: 445117.66335159 east: 480328.80428671 nsres: 30.02309444 ewres: 29.99245395 rows: 988 cols: 1174 cells: 1159912
Před dalším výpočtem je vhodné obrazové kanály družicové scény, které budou vstupovat do procesu klasifikace předzpracovat (aplikovat geometrické a atmosferické korekce). Minimálně nahradit místa bez informace hodnotou NULL - viz návod zde.
Dále je potřeba vytvořit obrazovou skupinu (i.group) slouženou ze vstupních kanálů - vynechán je panchromatický kanál 'B6'.
i.group group=B input=B1,B2,B3,B4,B5,B7
Ve wxGUI je pro vytvoření skupiny dostupný specializovaný interaktivní dialog (viz obr. níže).
Segmentace a statistika obrazových kanálů
Segmentaci obrazu provedeme příkazem i.segment.
i.segment group=B out=B_segment thresh=0.02 minsize=15
Vygenerované segmenty převedeme do vektorové reprezentace modulem r.to.vect a přidáme do atributové tabulky nové sloupce 'class' a 'brightness', zároveň odstraníme nepotřebný sloupec 'label'. Nepotřebné sloupečky odstraníme z atributové tabulky příkazem v.db.dropcolumn, nové potom v.db.addcolumn. Hodnoty nastavíme příkazem v.db.update.
# vektorizace segmentů
r.to.vect -v B_segment out=B_segment type=area
# modifikace atributové tabulky
v.db.dropcolumn B_segment col=label
v.db.addcolumn B_segment col="class int"
v.db.addcolumn B_segment col="brightness double precision"
v.db.update B_segment col=brightness value=0
Dále do atributové tabulky uložíme statistiku segmentů spojenou s jednotlivými kanály obrazové scény. Příklad pro první kanál:
v.rast.stats vector=B_segment raster=B1 column=B1
Jas ('brighness') vypočteme jako sumu středních hodnot jednotlivých kanálů:
v.db.update B_segment column=brightness qcolumn="brightness+B1_mean"
Po výpočtu můžeme nepotřebné sloupečky z atributové tabulky odstranit.
v.db.dropcolumn B_segment col=B1_n,B1_min,B1_max,B1_range,B1_cf_var,B1_sum
Statistiku můžeme spočítat dávkově pomocí grafického modeleru anebo pomocí jednoduchého skriptu v jazyce Python:
import grass.script as grass
maps = grass.read_command('i.group', group='B', flags='g', quiet=True).splitlines()
for map in maps:
band = map.split('@', 1)[0]
grass.message("Processing <%s>..." % band)
grass.run_command('v.rast.stats', vector='B_segment',
raster=band, column=band)
grass.run_command('v.db.update', map='B_segment',
column='brightness', qcolumn='brightness+%s_mean' % band)
grass.run_command('v.db.dropcolumn', map='B_segment',
column='{band}_n,{band}_min,{band}_max,' \
'{band}_range,{band}_cf_var,{band}_sum'.format(band=band))
Ve výsledku atributová tabulka obsahuje pro každý obrazový kanál informaci o střední hodnotě ('mean'), směrodatné odchylce ('stddev') a varianci nebo-li rozptylu ('variance') např.
B1_mean|72.6521739130435 B1_stddev|3.0447202434803 B1_variance|9.27032136105877 B2_mean|32.6086956521739 B2_stddev|1.99432085933559 B2_variance|3.97731568998105 ...
Analýza textury panchromatického kanálu
Analýzu textury provedeme příkazem r.texture:
r.texture input=B6 prefix=B6_text method=entr,idm
Jedné metodě odpovídá jedna rastrová mapa, v našem případě vzniknou dvě rastrové mapy: 'B6_text_IDM' a 'B6_text_Entr'.
Na základě těchto dvou rastrů aktualizujeme atributovou tabulku segmentů, zároveň odstraníme nepotřebné sloupečky. Zůstanou pouze: 'mean', 'stddev' a 'variance'. Statistiku do atributové tabulky uložíme příkazem v.rast.stats, nepotřebné sloupečky odstraníme v.db.dropcolumn.
# pro texturu IDM
v.rast.stats vector=B_segment raster=B6_text_IDM col=B6_text_IDM
v.db.dropcolumn B_segment col=B6_text_IDM_n,B6_text_IDM_min,B6_text_IDM_max,B6_text_IDM_range,B6_text_IDM_cf_var,B6_text_IDM_sum
# pro texturu Entr
v.rast.stats vector=B_segment raster=B6_text_Entr col=B6_text_Entr
v.db.dropcolumn B_segment col=B6_text_Entr_n,B6_text_Entr_min,B6_text_Entr_max,B6_text_Entr_range,B6_text_Entr_cf_var,B6_text_Entr_sum
Geometrické vlastnosti segmentů obrazu
V dalším kroku uložíme do atributové tabulky informace spojené s geometrií segmentů: jejich plochu ('area'), obvod ('perimeter'), kompaktnost ('compact') a fraktálovou dimenzi ('fd'). Sloupečky přidáme příkazem v.db.addcolumn anebo pomocí správce atributových dat, odvozené hodnoty z geometrie nahrajeme do atributové tabulky příkazem v.to.db.
# přidání nových sloupečků do atributové tabulky
v.db.addcolumn B_segment col="area int, perimeter int, compact double precision, fract_dim double precision" --q
# výpočet hodnot odvozených z geometrie segmentů
v.to.db B_segment op=area col=area
v.to.db B_segment op=perimeter col=perimeter
v.to.db B_segment op=compact col=compact
v.to.db B_segment op=fd col=fract_dim
Výsledek si ověříme (v.db.select nebo ve správci atributových dat) jednoduchým atributovým dotazem - vybereme segmenty s výměrou větší než 3 km2:
v.db.select B_segment col=cat,area,perimeter,compact,fract_dim where="area > 3e6"
cat|area|perimeter|compact|fract_dim 645628|3866602|35710|5.122953|1.382287 885913|3072390|37568|6.046235|1.410357
Volba trénovacích ploch
Volbu trénovacích ploch (resp. segmentů) provedeme v digitalizačním nástroji přiřazením příslušného kódu reprezentujího danou třídu (atribut 'class').
Trénovací segmenty vybereme a uložíme do nové vektorové mapy příkazem v.extract:
v.extract input=B_segment where="class > 0" output=tr_areas
Nastavíme vhodnou tabulku barev tak aby byly jednotlivé trénovací plochy dobře rozlišitelné. Interaktivní nástroj pro nastavení tabulky barev je k dispozici z menu wxGUI Vector → Manage colors → Manage color rules interactively.
Aplikace klasifikátoru
V našem případě použijeme klasifikátor DLDA implementovaný v rámci modulu v.class.mlpy.
Poznámka: Tento modul vyžaduje knihovnu mlpy.
v.class.mlpy input=B_segment training=tr_areas
Na závěř převedeme výsledek klasifikace do rastrové reprezentace pomocí v.to.rast:
v.to.rast input=B_segment output=B_obj_class use=attr attrcolumn=class
