153ZODH / 12. 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ů
Multitemporální analýza
K dispozici jsou následující data (lokace nc_spm_08):
- 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.
g.mapsets addmapset=modis2002lst
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.
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
# 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
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 sys
import time
import grass.script as grass
def check_monitor():
try:
monitor = grass.read_command('d.mon', flags = 'p').split(':')[1].strip()
except IndexError:
monitor = None
if not monitor:
grass.run_command('d.mon', start = 'x0')
else:
grass.run_command('d.mon', select = monitor)
def main():
if len(sys.argv) != 2:
print __doc__
sys.exit(1)
day = sys.argv[1]
mapset = 'modis2002lst'
# spustit graficky monitor
check_monitor()
# pridat mapset do vyhledavaci cestu
grass.run_command('g.mapsets',
addmapset = mapset)
for map in grass.mlist_grouped(type = 'rast',
mapset = mapset,
pattern = '*%s$' % day)[mapset]:
grass.message("Map: %s" % map)
grass.run_command('g.region',
rast = map)
grass.run_command('d.rast.leg',
quiet = True,
map = map)
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)
time.sleep(3)
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
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 9
"""
import sys
import time
import grass.script as grass
def check_monitor():
try:
monitor = grass.read_command('d.mon', flags = 'p').split(':')[1].strip()
except IndexError:
monitor = None
if not monitor:
grass.run_command('d.mon', start = 'x0')
else:
grass.run_command('d.mon', select = monitor)
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 = int(sys.argv[1])
mapset = 'modis2002lst'
# spustit graficky monitor
check_monitor()
# pridat mapset do vyhledavaci cestu
grass.run_command('g.mapsets',
addmapset = mapset)
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)
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)
time.sleep(1)
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
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
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
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()