153ZODH / 12. cvičení

← předchozí cvičení • stránky předmětu • další cvičení →

Multitemporální analýza

Osnova

Komplexní analýza obrazových dat je založena na monitorování zemského povrchu a atmosféry v různých prostorových rozlišeních a časových intervalech. Zaměříme se na analýzu časové řady dat MODIS (Moderate Resolution Imaging Spectroradiometer, satelit Terra a Aqua NASA).

Seznam příkazů

• g.mapsets
• g.list
• g.mlist
• d.rast.leg
• d.vect
• r.series
• r.mapcalc
• r.series
• r.colors

Multitemporální analýza

K dispozici jsou následující data (lokace nc_spm_08), ke stažení zde:

• Aqua (nosič, čas přeletu 1:30h)
• Terra (10:30h)
• Aqua (13:30h)
• Terra (22:30h)

Nejprve přidáme mapset obsahující data MODIS do vyhledávací cesty.

Z příkazové řádky

 g.mapsets addmapset=modis2002lst

či z wxGUI Settings → GRASS working environment → Mapset access

Nastavení vyhledávací cesty ve wxGUI

 g.list type=rast mapset=modis2002lst

V kombinaci se standardními Unixovými nástroji můžeme určit počet obrazových vrstev v mapsetu, seznam dnů pro které máme data k dispozici a pod.

Poznámka: Řetězení příkazů funguje pouze z příkazové řádky terminálu, nikoliv z wxGUI.

 g.mlist type=rast mapset=modis2002lst | wc -l

1008

 g.mlist type=rast mapset=modis2002lst | sed 's/[0-9]//g' | uniq

aqua_lst_day
aqua_lst_night
terra_lst_day
terra_lst_night

 g.mlist type=rast mapset=modis2002lst | sed -e 's/[a-z_]//g' | sort -n | uniq

20020101
20020102
...
20021230
20021231

Vybraná data zobrazíme.

 # noc, 22:30h, teplota ve stupních Celsia
 g.region rast=terra_lst_night20020921
 d.rast.leg map=terra_lst_night20020921
 d.vect map=us_natlas_hydrogp type=boundary color=blue
 d.vect map=us_natlas_urban type=boundary color=brown
 d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12

Obr. č.1: Teplota povrchu - 21.09.02, noc 22:30h, terra_lst_night20020921 (bílá místa představují odfiltrované mraky)

 # den, 10:30h
 d.rast.leg map=terra_lst_night20020922
 d.vect map=us_natlas_hydrogp type=boundary color=blue
 d.vect map=us_natlas_urban type=boundary color=brown
 d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12

Obr. č.2: Teplota povrhu - 22.09.02, noc 22:30h, terra_lst_night20020922

Příklad skriptu (Bash) pro vizualizaci dat za daný den.

#!/bin/sh

if [ -z "$1" ] ; then
    echo "Zadejte datum, napr. '$0 20020704'"
    exit 1
fi

day=$1

if [ `d.mon -p | cut -d':' -f2 | tr -d ' '` != 'x0' ] ; then
    d.mon x0 --q
fi

for map in `g.mlist type=rast mapset=modis2002lst pat="*${day}$"`; do
    g.message "Map: $map"
    g.region rast=$map
    d.rast.leg map=$map --q
    d.vect map=us_natlas_hydrogp type=boundary color=blue --q
    d.vect map=us_natlas_urban type=boundary color=brown --q
    d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12 --q

    sleep 3
done

exit 0

Verze skriptu pro Python:

#!/usr/bin/python
 
"""
Zadejte datum.
 
Napr.
./modis.py 20020704
"""
 
import os
import sys
import time
 
import grass.script as grass
 
def main():
    if len(sys.argv) != 2:
        print __doc__
        sys.exit(1)
 
    day = sys.argv[1]
    mapset = 'modis2002lst'
    monitor = 'cairo'
    
    # spustit graficky monitor
    grass.run_command('d.mon', start = monitor, width = 1024, height = 768)
 
    # vykreslit data
    for map in grass.mlist_grouped(type = 'rast',
                                   pattern = '*%s$' % day)[mapset]:
        map += '@' + mapset
        grass.message("Map: %s" % map)
        grass.run_command('g.region',
                          rast = map)
        grass.run_command('d.rast.leg',
                          quiet = True,
                          map = map)
        os.environ['GRASS_FRAME'] = "0.000000,768.000000,0.000000,716.800000"
        grass.run_command('d.vect',
                          quiet = True,
                          map = 'us_natlas_hydrogp',
                          type = 'boundary',
                          color = 'blue')
        grass.run_command('d.vect',
                          quiet = True,
                          map = 'us_natlas_urban',
                          type = 'boundary',
                          color = 'brown')
        grass.run_command('d.vect',
                          quiet = True,
                          map = 'us_natlas_urban',
                          display = 'attr',
                          attrcol = 'NAME',
                          type = 'centroid',
                          where = 'AREA >= 0.002',
                          lsize = 12)
        del os.environ['GRASS_FRAME']
 
        time.sleep(3)
 
    # vypnout graficky monitor
    grass.run_command('d.mon', stop = monitor)
    
    return 0
 
if __name__ == "__main__":
    sys.exit(main())

Rozdíl LST za 1 den vypočteme pomocí r.mapcalc.

 r.mapcalc "diff_den_noc = terra_lst_day20020922 - terra_lst_day20020921"
 d.rast.leg map=diff_den_noc
 d.vect map=us_natlas_hydrogp type=boundary color=blue
 d.vect map=us_natlas_urban type=boundary color=brown
 d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12

Obr. č.3: Rozdíl teplot povrchu za 1 den

Závěr: Na rozdíl od zastavěných ploch se teplotně stabilní jeví vodní plochy.

Indikátory časových řad

Data jsou již předzpracována - digitální hodnoty odpovídající teplotě mraků byly odfiltrovány. Během výpočtu se hodnoty NULL mohou výrazněji promítnout do výsledku, např. při určování průměrné teploty povrchu.

Vypočteme denní LST průměry pro září 2002, modul r.series.

for d in `seq 1 30` ; do
    DAY=`echo $d | awk '{printf "%02d\n", $1}'`
    LIST=`g.mlist type=rast pattern="*lst*200209${DAY}" separator=","`
    g.message "Den: $DAY ($LIST)"
    g.region rast=$LIST
    r.series -n input=$LIST output=lst_200209${DAY}_avg method=average --o --q
    d.rast.leg map=lst_200209${DAY}_avg --q
    d.vect map=us_natlas_hydrogp type=boundary color=blue --q
    d.vect map=us_natlas_urban type=boundary color=brown --q
    d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12 --q
    
    sleep 1
done

Jako ukázka korespondující skript v Pythonu:

 #!/usr/bin/python
 
"""
Zadejte mesic (1-12)
 
Napr.
./modis.py 09
"""

import os
import sys
import time
 
import grass.script as grass
 
def get_num_days(month):
    if month in (1, 3, 5, 7, 8, 10, 12):
        return 32
    elif month == 2:
        return 29
 
    return 31
 
def main():
    if len(sys.argv) != 2:
        print __doc__
        sys.exit(1)
 
    month = sys.argv[1]
    mapset = 'modis2002lst'
    monitor = 'cairo'
    
    # spustit graficky monitor
    grass.run_command('d.mon', start = monitor, width = 1024, height = 768)
 
    # pridat mapset do vyhledavaci cestu
    grass.run_command('g.mapsets', mapset = mapset, operation = 'add', quiet = True)
  
    for day in range(1, get_num_days(month)):
        day = '%02d' % day
        maps = ','.join(grass.mlist_grouped(type = 'rast',
                                            pattern = "*lst*2002%s%s" % (month, day))[mapset])
        omap = 'lst_2002%s%s_avg' % (month, day)
 
        grass.message('Den: %s (%s)' % (day, maps))
        grass.run_command('g.region',
                          rast = maps)
        grass.run_command('r.series',
                          quiet = True,
                          overwrite = True,
                          flags = 'n',
                          input = maps,
                          output = omap,
                          method = 'average')
        grass.run_command('d.rast.leg',
                          quiet = True,
                          map = omap)
        os.environ['GRASS_FRAME'] = "0.000000,768.000000,0.000000,716.800000"
        grass.run_command('d.vect',
                          quiet = True,
                          map = 'us_natlas_hydrogp',
                          type = 'boundary',
                          color = 'blue')
        grass.run_command('d.vect',
                          quiet = True,
                          map = 'us_natlas_urban',
                          type = 'boundary',
                          color = 'brown')
        grass.run_command('d.vect',
                          quiet = True,
                          map = 'us_natlas_urban',
                          display = 'attr',
                          attrcol = 'NAME',
                          type = 'centroid',
                          where = 'AREA >= 0.002',
                          lsize = 12)
        del os.environ['GRASS_FRAME']

        time.sleep(3)
 
    # vypnout graficky monitor
    grass.run_command('d.mon', stop = monitor)

    return 0
 
if __name__ == "__main__":
    sys.exit(main())

V důsledku NULL hodnot ve vstupních datech některé výstupní vrstvy obsahují částečně nebo zcela hodnoty NULL. Odstranění tohoto efektu by vyžadovalo komplexní rekonstrukci časové řady.

Podobně vytvoříme vrstvu průměrné LST za měsíc září 2002.

 LIST=`g.mlist type=rast pattern="*lst*200209[0-9][0-9]$" sep=","`
 g.region rast=$LIST
 r.series input=$LIST output=lst_200209_avg method=average
 r.colors map=lst_200209_avg color=gyr

Obr č.4: Průměrná teplota povrchu září 2002

Určíme počet validních pixelů (not NULL) v časové řadě.

 r.series input=$LIST output=lst_200209_avg_count method=count
 r.colors map=lst_200209_avg_count color=gyr

 d.rast.leg map=lst_200209_avg_count
 d.vect map=us_natlas_hydrogp type=boundary color=blue
 d.vect map=us_natlas_urban type=boundary color=brown
 d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12

Obr č.5: Počet validních pixelů (not NULL)

Na základě počtu validních pixelů nastavíme filtr.

 r.mapcalc "lst_200209_avg_filt = if(lst_200209_avg_count > 5, lst_200209_avg, null())"

 d.rast.leg map=lst_200209_avg_filt
 d.vect map=us_natlas_hydrogp type=boundary color=blue
 d.vect map=us_natlas_urban type=boundary color=brown
 d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12

Obr č.6: Průměrná teplota povrchu září 2002 (odfiltrovány místa s nízkým počtem validních pixelů

Zastavěné plochy přirozeně vykazují vyšší průměrnou teplotu v porovnání např. s vodními plochami.

Ukázka skriptu Python

#!/usr/bin/python

"""
Zadejte mesic (1-12)

Napr.
./modis.py 9
"""

import sys
import time

import grass.script as grass

def main():
    mapset = 'modis2002lst'
    
    mapsets = grass.parse_command('g.mapsets',
                                  flags = 'p',
                                  fs = '\n').keys()
    
    if mapset not in mapsets:
        # pridat mapset do vyhledavaci cesty
        grass.run_command('g.mapsets',
                          addmapset = mapset)

    maps_avg = list()
    for month in range(1, 13):
        month = '%02d' % month
        maps_avg.append('lst_2002%s_avg' % month)
        
        grass.message('Processing month %s...' % month)
        
        maps = ','.join(grass.mlist_grouped(type = 'rast',
                                            mapset = mapset,
                                            pattern = '*lst*2002%s[0-9][0-9]$' % month)[mapset])
        
        grass.run_command('g.region',
                          rast = maps)
        grass.run_command('r.series',
                          quiet = True,
                          overwrite = True,
                          input = maps,
                          output = maps_avg[-1],
                          method = 'average')
        
        grass.run_command('r.colors',
                          quiet = True,
                          map = maps_avg[-1],
                          color = 'gyr')
        stats = grass.read_command('r.univar',
                                quiet = True,
                                flags = 'g',
                                map = maps_avg[-1])
        stats = grass.parse_key_val(stats)
        for key in stats.iterkeys():
            if key in ('min', 'max', 'mean'):
                value = int(round(float(stats[key])))
                print '%s = %02dC' % (key, value)
        
    grass.run_command('g.mremove',
                      rast = ','.join(maps_avg))
    
    if mapsets not in mapsets:
        # odebrat mapset z vyhledavaci cesty
        grass.run_command('g.mapsets',
                          removemapset = mapset)
    
if __name__ == "__main__":
    main()