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).

Skripty v jazyku Python používají knihovnu pyGRASS.

Seznam příkazů

• g.mapsets
• g.list
• g.mlist
• d.rast.leg
• d.vect
• r.series
• r.mapcalc
• r.series
• r.colors
• r.timestamp
• t.create
• t.register
• t.info

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 mapset=modis2002lst operation=add

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

Nastavení vyhledávací cesty ve wxGUI

 g.list type=rast mapset=modis2002lst

V jazyku Python:

from grass.pygrass.gis import Mapset

print Mapset('modis2002lst').glist('rast')

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

V jazyku Python:

print len(Mapset('modis2002lst').glist('rast'))

1008

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

V jazyku Python:

import re

print set(map(lambda x: re.sub(r'[0-9]', '', x), Mapset('modis2002lst').glist('rast')))

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

V jazyku Python:

for date in sorted(set(map(lambda x: re.sub(r'[a-z_]', '', x), Mapset('modis2002lst').glist('rast')))):
    print date

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

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

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

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

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

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

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
 
from subprocess import PIPE

from grass.pygrass.gis import Location, Mapset
from grass.pygrass.modules import Module
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r

from grass.script import parse_key_val
 
def main():
    modis_mapset = 'modis2002lst'
 
    if modis_mapset not in Location().mapsets():
        # pridat mapset do vyhledavaci cesty
        g.mapsets(mapset = modis_mapset, operation = 'add')
    
    maps_avg = list()
    for month in range(1, 13):
        maps_avg.append('lst_2002%02d_avg' % month)
        
        g.message('Processing month %02d...' % month)
        maps = ','.join(sorted(Mapset(modis_mapset).glist('rast',
                                                          pattern = '*lst*2002%02d[0-9][0-9]' % month)))

        g.region(rast = maps)
        r.series(quiet = True,
                 overwrite = True,
                 input = maps,
                 output = maps_avg[-1],
                 method = 'average')
 
        r.colors(quiet = True,
                 map = maps_avg[-1],
                 color = 'gyr')

        univar = Module("r.univar")
        univar(quiet = True,
               flags = 'g',
               map = maps_avg[-1], stdout_=PIPE)
        stats = parse_key_val(univar.outputs.stdout)
        
        for key in stats.iterkeys():
            if key in ('min', 'max', 'mean'):
                value = int(round(float(stats[key])))
                print '%s = %02dC' % (key, value)
    
    g.mremove(rast = ','.join(maps_avg), flags = 'f', quiet=True)
 
if __name__ == "__main__":
    main()

Výstup:

Processing month 01...
min = 06C
max = 14C
mean = 10C
Processing month 02...
min = 07C
max = 14C
mean = 11C
Processing month 03...
min = 09C
max = 29C
mean = 18C
...

Vytvoření a analýza časoprostorových dat

Tato část je založena na Temporal framework v systému GRASS.

Nejprve vytvoříme databázi pro registraci časoprostorového datasetu.

 t.create output=modis title="MODIS 2002" desc="Ukazkovy casoprostorovy dataset MODIS"

V dalším kroku zaregistrujeme MODIS data. Tato data nelze registrovat automaticky neboť u nich chybí časová značka. Např.

r.info map=terra_lst_night20021231

...
 | Timestamp: none                                                            
...

Chybějící časovou značku můžeme dopnit pomocí následujícího skriptu:

#!/usr/bin/env python

#%module
#% description: Set timestamps for MODIS LTS raster maps
#%end

#%option
#% key: from
#% description: From date (YYYY-MM-DD)
#% required: yes
#%end

#%option
#% key: to
#% description: To date (YYYY-MM-DD)
#% required: yes
#%end

#%option G_OPT_F_OUTPUT
#% description: Path to the map file
#% required: no
#%end

import os
import sys
import datetime
import calendar

from grass.script.core import parser, fatal, message, warning
from grass.pygrass.gis import Mapset
from grass.pygrass.modules.shortcuts import raster as r

def get_date(date):
    try:
        y, m, d = map(int, date.split('-'))
    except:
        fatal("Invalid date format '%s'" % date)
    
    return datetime.date(y, m, d)

def get_month(date):
    return calendar.month_abbr[date.month]

def main():
    fd = None
    if options['output']:
        fd = open(options['output'], 'w')
    
    date = from_date = get_date(options['from'])
    to_date = get_date(options['to'])
    while date < to_date:
        message("%s..." % str(date))
        
        map_timestamp = '%s %s %s' % (date.day, get_month(date), date.year)
        map_date = str(date).replace('-', '')
        
        maps = Mapset().glist('rast', pattern='*%s' % map_date)
        if len(maps) != 4:
            warning("Not complete (%d)" % len(maps))
        
        for rast, time in [("aqua_lst_night", "1:30"),
                           ("terra_lst_day", "10:30"),
                           ("aqua_lst_day", "13:30"),
                           ("terra_lst_night", "22:30")]:
            map_name = rast + map_date
            if map_name not in maps:
                continue
            
            r.timestamp(map = map_name, date = map_timestamp + ' ' + time)
            if fd:
                fd.write(map_name + os.linesep)
            
        date += datetime.timedelta(days=1)
    
    if fd:
        fd.close()
        
if __name__ == "__main__":
    options, flags = parser()
    sys.exit(main())

r.stampmodis.py from=2002-01-01 to=2003-01-01 out=maps.txt

r.info map=terra_lst_night20021231

...
 | Timestamp: 31 Dec 2002 22:30                                                             
...

Poznámka: Skript vygeneruje soubor se seznamem map, který použijeme jako vstup do modulu t.register.

Po této úpravě můžeme data úspěšně zaregistrovat

 t.register input=modis file=maps.txt

Základní metadata časoprostorového datasetu vypíšeme pomocí modulu t.info:

t.info modis

 +-------------------- Space Time Raster Dataset -----------------------------+
 |                                                                            |
 +-------------------- Basic information -------------------------------------+
 ...
 | Temporal type: ............. absolute
 ...
 | Semantic type:.............. mean
 +-------------------- Absolute time -----------------------------------------+
 | Start time:................. 2002-01-01 10:30:00
 | End time:................... 2002-12-31 22:30:00
 | Granularity:................ 180 minutes
 | Temporal type of maps:...... point
 +-------------------- Spatial extent ----------------------------------------+
 | North:...................... 279073.975466
 | South:...................... 113673.975466
 | East:.. .................... 798143.311797
 | West:....................... 595143.311797
 | Top:........................ 0.0
 | Bottom:..................... 0.0
 +-------------------- Metadata information ----------------------------------+
 | Raster register table:...... modis_modis2002lst_1_raster_register
 | North-South resolution min:. 200.0
 | North-South resolution max:. 200.0
 | East-west resolution min:... 200.0
 | East-west resolution max:... 200.0
 | Minimum value min:.......... -1.49
 | Minimum value max:.......... 30.63
 | Maximum value min:.......... -0.29
 | Maximum value max:.......... 48.23
 | Aggregation type:........... None
 | Number of registered maps:.. 1008
 |
 | Title:
 | MODIS 2002
 | Description:
 | Ukazkovy casoprostorovy dataset MODIS
 | Command history:
 | # 2014-04-22 21:13:08 
 | t.create output="modis" title="MODIS 2002"
 |     desc="Ukazkovy casoprostorovy dataset MODIS"
 | # 2014-04-22 21:18:05 
 | t.register input="modis" file="maps.txt"
 | 
 +----------------------------------------------------------------------------+

Kontrola časové topologie

t.topology modis

 +-------------------- Basic information -------------------------------------+
 | Id: ........................ modis@modis2002lst_1
 | Name: ...................... modis
 | Mapset: .................... modis2002lst_1
 | Creator: ................... landa
 | Temporal type: ............. absolute
 | Creation time: ............. 2014-04-22 21:13:08.494237
 | Modification time:.......... 2014-04-22 21:18:05.749742
 | Semantic type:.............. mean
 +-------------------- Temporal topology -------------------------------------+
 | Is subset of dataset: ...... False
 | Temporal topology is: ...... valid
 | Number of intervals: ....... 0
 | Invalid time stamps: ....... 0
 | Number of points: .......... 1008
 | Number of gaps: ............ 1007
 | Granularity: ............... 180 minutes
 +-------------------- Topological relations ---------------------------------+
 | Overlaps: .................. 0
 | Overlapped: ................ 0
 | Finishes: .................. 0
 | Started: ................... 0
 | Follows: ................... 0
 | Contains: .................. 0
 | Equal:...................... 0
 | Finished: .................. 0
 | Precedes: .................. 0
 | Starts: .................... 0
 | During: .................... 0
 +----------------------------------------------------------------------------+

Základní dotazy

Modul t.rast.list umožnuje vypsat z datasetu rastrové mapy odpovídající dané podmínkce, např. všechny snímky z měsíce března:

t.rast.list input=modis order=start_time where="start_time > '2002-03-01' and start_time < '2002-04-01'"

terra_lst_day20020312	modis2002lst_1	2002-03-12 10:30:00	None
terra_lst_night20020312	modis2002lst_1	2002-03-12 22:30:00	None
terra_lst_day20020314	modis2002lst_1	2002-03-14 10:30:00	None
...
terra_lst_night20020330	modis2002lst_1	2002-03-30 22:30:00	None
terra_lst_day20020331	modis2002lst_1	2002-03-31 10:30:00	None
terra_lst_night20020331	modis2002lst_1	2002-03-31 22:30:00	None

Základní statistiku rastrových map poskutuje modul t.rast.univar:

t.rast.univar input=modis where="start_time > '2002-03-01' and start_time < '2002-04-01'"

terra_lst_night20020312@modis2002lst_1|2002-03-12 22:30:00|None|7.12178469216311|6.61000000000001|11.45|...
...
terra_lst_day20020331@modis2002lst_1|2002-03-31 10:30:00|None|21.9590689981096|20.45|27.51|...

Vizualizace časoprostorových dat

g.gui.animation strds=modis