153ZODH / 12. cvičení: Porovnání verzí
m (155zddp) |
|||
(Není zobrazeno 152 mezilehlých verzí od stejného uživatele.) | |||
Řádek 1: | Řádek 1: | ||
{{ | {{Zastaralé|155ZDDP}} | ||
{{upravit}} | |||
{{Cvičení|153ZODH|12|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 {{wikipedia|MODIS|lang=en}} (''Moderate Resolution Imaging Spectroradiometer'', satelit Terra a Aqua NASA). | |||
Skripty v jazyku [[Python]] používají knihovnu [http://grasswiki.osgeo.org/wiki/Python/pygrass pyGRASS]. | |||
=== Seznam příkazů === | |||
* {{GrassPrikaz|g.mapsets}} | |||
* {{GrassPrikaz|g.list}} | |||
* {{GrassPrikaz|g.mlist}} | |||
* {{GrassPrikaz|d.rast.leg}} | |||
* {{GrassPrikaz|d.vect}} | |||
* {{GrassPrikaz|r.series}} | |||
* {{GrassPrikaz|r.mapcalc}} | |||
* {{GrassPrikaz|r.series}} | |||
* {{GrassPrikaz|r.colors}} | |||
* {{GrassPrikaz|r.timestamp}} | |||
* {{GrassPrikaz|t.create}} | |||
* {{GrassPrikaz|t.register}} | |||
* {{GrassPrikaz|t.info}} | |||
== Multitemporální analýza == | |||
K dispozici jsou následující data (lokace [http://grass.osgeo.org/sampledata/north_carolina/nc_spm_08_grass7.tar.gz nc_spm_08]), ke stažení [http://grass.osgeo.org/sampledata/north_carolina/nc_spm_modis2002lst_vector.tar.gz zde]: | |||
* Aqua (1:30h) | * Aqua (nosič, čas přeletu 1:30h) | ||
* Terra (10:30h) | * Terra (10:30h) | ||
* Aqua (13:30h) | * Aqua (13:30h) | ||
* Terra (22:30h) | * Terra (22:30h) | ||
=== | Nejprve přidáme mapset obsahující data MODIS do vyhledávací cesty. | ||
Z příkazové řádky | |||
<source lang="bash"> | |||
g.mapsets mapset=modis2002lst operation=add | |||
</source> | |||
či z wxGUI {{menu|Settings|GRASS working environment|Mapset access}} | |||
{{fig|wxgui-mapsets|Nastavení vyhledávací cesty ve wxGUI|size=300}} | |||
<source lang="bash"> | |||
g.list type=rast mapset=modis2002lst | |||
</source> | |||
V jazyku Python: | |||
<source lang=python> | |||
from grass.pygrass.gis import Mapset | |||
print Mapset('modis2002lst').glist('rast') | |||
</source> | |||
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 {{grassPrikaz|wxGUI}}. | |||
<source lang="bash"> | |||
g.mlist type=rast mapset=modis2002lst | wc -l | |||
</source> | |||
V jazyku Python: | |||
<source lang=python> | |||
print len(Mapset('modis2002lst').glist('rast')) | |||
</source> | |||
1008 | |||
<source lang="bash"> | |||
g.mlist type=rast mapset=modis2002lst | sed 's/[0-9]//g' | uniq | |||
</source> | |||
V jazyku Python: | |||
<source lang=python> | |||
import re | |||
print set(map(lambda x: re.sub(r'[0-9]', '', x), Mapset('modis2002lst').glist('rast'))) | |||
</source> | |||
aqua_lst_day | |||
aqua_lst_night | |||
terra_lst_day | |||
terra_lst_night | |||
<source lang="bash"> | |||
g.mlist type=rast mapset=modis2002lst | sed -e 's/[a-z_]//g' | sort -n | uniq | |||
</source> | |||
V jazyku Python: | |||
<source lang=python> | |||
for date in sorted(set(map(lambda x: re.sub(r'[a-z_]', '', x), Mapset('modis2002lst').glist('rast')))): | |||
print date | |||
</source> | |||
20020101 | |||
20020102 | |||
... | |||
20021230 | |||
20021231 | |||
Vybraná data zobrazíme. | |||
<source lang="bash"> | |||
# noc, 22:30h, teplota ve stupních Celsia | # noc, 22:30h, teplota ve stupních Celsia | ||
d.rast.leg map= | 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_hydrogp type=boundary color=blue | ||
d.vect map=us_natlas_urban type=boundary | d.vect map=us_natlas_urban type=boundary color=brown | ||
d.vect map=us_natlas_urban | d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12 | ||
</source> | |||
{{fig|ZOD-cv12-terra-first|Teplota povrchu - 21.09.02, noc 22:30h, terra_lst_night20020921 (bílá místa představují odfiltrované mraky)}} | |||
<source lang="bash"> | |||
# den, 10:30h | # den, 10:30h | ||
d.rast.leg map= | 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 | |||
</source> | |||
{{fig|ZOD-cv12-terra-second|Teplota povrhu - 22.09.02, noc 22:30h, terra_lst_night20020922}} | |||
Příklad skriptu (Bash) pro vizualizaci dat za daný den. | |||
<source lang="bash"> | |||
#!/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 | |||
</source> | |||
Verze skriptu pro Python: | |||
<source lang="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()) | |||
</source> | |||
Rozdíl LST za 1 den vypočteme pomocí {{GrassPrikaz|r.mapcalc}}. | |||
<source lang="bash"> | |||
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_hydrogp type=boundary color=blue | ||
d.vect map=us_natlas_urban type=boundary | d.vect map=us_natlas_urban type=boundary color=brown | ||
d.vect map=us_natlas_urban | d.vect map=us_natlas_urban display=attr attrcol=NAME type=centroid where="AREA >= 0.002" lsize=12 | ||
</source> | |||
{{fig|ZOD-cv12-terra-diff|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 {{GrassPrikaz|r.series}}. | |||
<source lang="bash"> | |||
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 | |||
</source> | |||
Jako ukázka korespondující skript v Pythonu: | |||
<source lang="python"> | |||
#!/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()) | |||
</source> | |||
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. | |||
<source lang="bash"> | |||
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 | |||
</source> | |||
{{fig|ZOD-cv12-prumerna-teplota-zari|Průměrná teplota povrchu září 2002}} | |||
Určíme počet validních pixelů (not NULL) v časové řadě. | |||
<source lang="bash"> | |||
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 | |||
</source> | |||
{{fig|ZOD-cv12-prumerna-teplota-zari-pocet|Počet validních pixelů (not NULL)}} | |||
Na základě počtu validních pixelů nastavíme filtr. | |||
<source lang="bash"> | |||
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 | |||
</source> | |||
{{fig|ZOD-cv12-prumerna-teplota-zari-filtr|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 == | |||
<source lang="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() | |||
</source> | |||
Výstup: | |||
<pre> | |||
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 | |||
... | |||
</pre> | |||
== Vytvoření a analýza časoprostorových dat == | |||
Tato část je založena na {{grassPrikaz|temporalintro|Temporal framework}} v systému GRASS. | |||
Nejprve vytvoříme databázi pro registraci časoprostorového datasetu. | |||
<source lang=bash> | |||
t.create output=modis title="MODIS 2002" desc="Ukazkovy casoprostorovy dataset MODIS" | |||
</source> | |||
V dalším kroku zaregistrujeme MODIS data. Tato data nelze registrovat automaticky neboť u nich chybí časová značka. Např. | |||
<source lang=bash> | |||
r.info map=terra_lst_night20021231 | |||
</source> | |||
<pre> | |||
... | |||
| Timestamp: none | |||
... | |||
</pre> | |||
Chybějící časovou značku můžeme dopnit pomocí následujícího skriptu: | |||
<source lang=python line> | |||
#!/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()) | |||
</source> | |||
<source lang=bash> | |||
r.stampmodis.py from=2002-01-01 to=2003-01-01 out=maps.txt | |||
</source> | |||
<source lang=bash> | |||
r.info map=terra_lst_night20021231 | |||
</source> | |||
<pre> | |||
... | |||
| Timestamp: 31 Dec 2002 22:30 | |||
... | |||
</pre> | |||
''Poznámka:'' Skript vygeneruje soubor se seznamem map, který použijeme jako vstup do modulu {{grassPrikaz|t.register}}. | |||
Po této úpravě můžeme data úspěšně zaregistrovat | |||
<source lang=bash> | |||
t.register input=modis file=maps.txt | |||
</source> | |||
Základní metadata časoprostorového datasetu vypíšeme pomocí modulu {{grassPrikaz|t.info}}: | |||
<source lang=python> | |||
t.info modis | |||
</source> | |||
<pre> | |||
+-------------------- 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" | |||
| | |||
+----------------------------------------------------------------------------+ | |||
</pre> | |||
==== Kontrola časové topologie ==== | |||
Souhrné informace o časové topologii vystiskne příkaz {{grassPrikaz|t.topology}}: | |||
<source lang=python> | |||
t.topology modis | |||
</source> | |||
<pre> | |||
+-------------------- 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 | |||
... | |||
</pre> | |||
Pro podrobné informace je nutné přidat přepínač <tt>-m</tt>: | |||
<source lang=bash> | |||
t.topology -m modis where="start_time >= '2002-07-15' and start_time < '2002-07-17'" | |||
</source> | |||
<pre> | |||
+-------------------- Raster Dataset ----------------------------------------+ | |||
| | | |||
+-------------------- Basic information -------------------------------------+ | |||
| Id: ........................ aqua_lst_night20020715@modis2002lst_yfsg | |||
| Name: ...................... aqua_lst_night20020715 | |||
| Mapset: .................... modis2002lst_yfsg | |||
| Creator: ................... martin | |||
| Temporal type: ............. absolute | |||
| Creation time: ............. 2014-04-29 18:34:29.485248 | |||
+-------------------- Absolute time -----------------------------------------+ | |||
| Start time:................. 2002-07-15 01:30:00 | |||
| End time:................... None | |||
+-------------------- Temporal Topology -------------------------------------+ | |||
| Next: ...................... terra_lst_day20020715@modis2002lst_yfsg | |||
+-------------------- Spatial extent ----------------------------------------+ | |||
| North:...................... 279073.975466 | |||
| South:...................... 113673.975466 | |||
| East:.. .................... 798143.311797 | |||
| West:....................... 595143.311797 | |||
| Top:........................ 0.0 | |||
| Bottom:..................... 0.0 | |||
+-------------------- Metadata information ----------------------------------+ | |||
| Datatype:................... DCELL | |||
| Number of columns:.......... 827 | |||
| Number of rows:............. 1015 | |||
| Number of cells:............ 839405 | |||
| North-South resolution:..... 200.0 | |||
| East-west resolution:....... 200.0 | |||
| Minimum value:.............. 12.85 | |||
| Maximum value:.............. 19.53 | |||
| Registered datasets ........ modis@modis2002lst_yfsg | |||
+----------------------------------------------------------------------------+ | |||
... | |||
+-------------------- Raster Dataset ----------------------------------------+ | |||
| | | |||
+-------------------- Basic information -------------------------------------+ | |||
| Id: ........................ terra_lst_night20020716@modis2002lst_yfsg | |||
| Name: ...................... terra_lst_night20020716 | |||
| Mapset: .................... modis2002lst_yfsg | |||
| Creator: ................... martin | |||
| Temporal type: ............. absolute | |||
| Creation time: ............. 2014-04-29 18:34:29.603983 | |||
+-------------------- Absolute time -----------------------------------------+ | |||
| Start time:................. 2002-07-16 22:30:00 | |||
| End time:................... None | |||
+-------------------- Temporal Topology -------------------------------------+ | |||
| Previous: .................. aqua_lst_day20020716@modis2002lst_yfsg | |||
+-------------------- Spatial extent ----------------------------------------+ | |||
| North:...................... 279073.975466 | |||
| South:...................... 113673.975466 | |||
| East:.. .................... 798143.311797 | |||
| West:....................... 595143.311797 | |||
| Top:........................ 0.0 | |||
| Bottom:..................... 0.0 | |||
+-------------------- Metadata information ----------------------------------+ | |||
| Datatype:................... DCELL | |||
| Number of columns:.......... 827 | |||
| Number of rows:............. 1015 | |||
| Number of cells:............ 839405 | |||
| North-South resolution:..... 200.0 | |||
| East-west resolution:....... 200.0 | |||
| Minimum value:.............. 18.95 | |||
| Maximum value:.............. 24.33 | |||
| Registered datasets ........ modis@modis2002lst_yfsg | |||
+----------------------------------------------------------------------------+ | |||
</pre> | |||
==== Základní dotazy ==== | |||
Modul {{grassPrikaz|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: | |||
<source lang=bash> | |||
t.rast.list input=modis order=start_time where="start_time > '2002-03-01' and start_time < '2002-04-01'" | |||
</source> | |||
<pre> | |||
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 | |||
</pre> | |||
Základní statistiku rastrových map poskutuje modul {{grassPrikaz|t.rast.univar}}: | |||
<source lang=bash> | |||
t.rast.univar input=modis where="start_time > '2002-03-01' and start_time < '2002-04-01'" | |||
</source> | |||
<pre> | |||
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|... | |||
</pre> | |||
==== Agregace ==== | |||
Určení statististiky teplot pro jednotlivé měsíce: | |||
<source lang=python> | |||
t.rast.aggregate input=modis out=modis_m basename=ag gra="1 months" | |||
</source> | |||
Vytvoří se dvanáct rastrových map v měsíční periodě: | |||
<source lang=python> | |||
t.rast.list modis_m order=start_time | |||
</source> | |||
<pre> | |||
ag_0 modis2002lst_1 2002-01-01 00:00:00 2002-02-01 00:00:00 | |||
ag_1 modis2002lst_1 2002-02-01 00:00:00 2002-03-01 00:00:00 | |||
ag_2 modis2002lst_1 2002-03-01 00:00:00 2002-04-01 00:00:00 | |||
ag_3 modis2002lst_1 2002-04-01 00:00:00 2002-05-01 00:00:00 | |||
ag_4 modis2002lst_1 2002-05-01 00:00:00 2002-06-01 00:00:00 | |||
ag_5 modis2002lst_1 2002-06-01 00:00:00 2002-07-01 00:00:00 | |||
ag_6 modis2002lst_1 2002-07-01 00:00:00 2002-08-01 00:00:00 | |||
ag_7 modis2002lst_1 2002-08-01 00:00:00 2002-09-01 00:00:00 | |||
ag_8 modis2002lst_1 2002-09-01 00:00:00 2002-10-01 00:00:00 | |||
ag_9 modis2002lst_1 2002-10-01 00:00:00 2002-11-01 00:00:00 | |||
ag_10 modis2002lst_1 2002-11-01 00:00:00 2002-12-01 00:00:00 | |||
ag_11 modis2002lst_1 2002-12-01 00:00:00 2003-01-01 00:00:00 | |||
</pre> | |||
Statistika pro všechny měsíce: | |||
<source lang=python> | |||
t.rast.univar modis_m | |||
</source> | |||
<pre> | |||
ag_0@modis2002lst_1|2002-01-01 00:00:00|2002-02-01 00:00:00|10.1761589696136|6.11111111111114|14.3009090909091|... | |||
ag_11@modis2002lst_1|2002-12-01 00:00:00|2003-01-01 00:00:00|7.50005058325719|4.3941666666667|10.3609523809524|... | |||
</pre> | |||
Statistika pouze pro červenec a srpen: | |||
<source lang=python> | |||
t.rast.univar modis_m where="start_time > '2002-07-01' and start_time < '2002-09-01'" | |||
</source> | |||
<pre> | |||
ag_6@modis2002lst_1|2002-07-01 00:00:00|2002-08-01 00:00:00|26.3466550685411|21.7652941176471|30.594|... | |||
ag_7@modis2002lst_1|2002-08-01 00:00:00|2002-09-01 00:00:00|26.7082631007969|22.2151851851852|31.4769565217392|... | |||
</pre> | |||
==== Výběr dat z časoprostorového datasetu ==== | |||
Vytvořit na základě výběru nový časoprostorový dataset umožňuje příkaz {{grassPrikaz|t.rast.extract}}. | |||
<source lang=bash> | |||
t.rast.extract input=modis where="start_time > '2002-03-01' and start_time < '2002-06-01'" output=modis_spring | |||
t.rast.extract input=modis where="start_time > '2002-06-01' and start_time < '2002-09-01'" output=modis_summer | |||
t.rast.extract input=modis where="start_time > '2002-09-01' and start_time < '2002-12-01'" output=modis_autumn | |||
t.rast.extract input=modis where="start_time > '2002-12-01' or start_time < '2002-03-01'" output=modis_winter | |||
</source> | |||
Vizualizace časové řady pomocí nástroje {{grassPrikaz|g.gui.timeline}}. | |||
<source lang=bash> | |||
g.gui.timeline inputs=modis_spring,modis_summer,modis_autumn,modis_winter | |||
</source> | |||
{{fig|g-gui-timeline-modis-4|Vizualizace čtyř časoprostorových datasetů na základě ročního období|size=640}} | |||
V následujících příkazech budeme sledovat trend změny teploty v jednotlivých ročních obdobích. K tomu použijeme modul {{grassPrikaz|t.rast.series}}. | |||
<source lang=python> | |||
t.rast.series input=modis_spring output=modis_spring_avg method=average | |||
t.rast.series input=modis_summer output=modis_summer_avg method=average | |||
t.rast.series input=modis_autumn output=modis_autumn_avg method=average | |||
t.rast.series input=modis_winter output=modis_winter_avg method=average | |||
</source> | |||
Vzniknou čtyři rastrové mapy zobrazující průměrné roční teploty v ročních obdobích. Průměrnou teplotu zjistíme pomocí modulu {{grassPrikaz|r.univar}}, příklad pro jaro: | |||
<source lang=python> | |||
r.univar modis_spring_avg | |||
... | |||
mean: 21.6125 | |||
... | |||
</source> | |||
{{fig|modis-spring-avg|Průměrné teploty pro roční období 'jaro' (tabulka barev 'celsius')}} | |||
==== Vizualizace časoprostorových dat ==== | |||
Vizualizovat data časoprostorových datasetů umožňuje {{grassPrikaz|g.gui.animation|wxGUI animační nástroj}}. | |||
<source lang=bash> | |||
g.gui.animation strds=modis | |||
</source> | |||
{{fig|g-gui-animation|wxGUI Animation Tool|size=640}} | |||
Mezi další užitečné nástroje patří {{grassPrikaz|g.gui.mapswipe|wxGUI Map Swipe}} | |||
<source lang=bash> | |||
t.rast.list modis_m where="start_time < '2002-03-01'" | |||
</source> | |||
<pre> | |||
ag_0 modis2002lst_yfsg 2002-01-01 00:00:00 2002-02-01 00:00:00 | |||
ag_1 modis2002lst_yfsg 2002-02-01 00:00:00 2002-03-01 00:00:00 | |||
</pre> | |||
<source lang=bash> | |||
g.gui.mapswipe first=ag_0 second=ag_1 | |||
</source> | |||
{{fig|g-gui-mapswipe-modis|Vizualizace agregovaných LTS dat pro první dva měsíce roku 2002|size=640}} | |||
a {{grassPrikaz|g.gui.timeline}}: | |||
--- | {{fig|g-gui-timeline-modis|Vizualizace časové osy agregovaných LTS dat|size=640}} | ||
{{ZOD}} | {{ZOD}} | ||
{{GRASS}} | |||
{{Python}} |
Aktuální verze z 5. 9. 2014, 18:49
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
![](/images/0/06/Wxgui-mapsets.png)
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
![](/images/e/e2/ZOD-cv12-terra-first.png)
# 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
![](/images/0/00/ZOD-cv12-terra-second.png)
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
![](/images/5/57/ZOD-cv12-terra-diff.png)
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
![](/images/a/a4/ZOD-cv12-prumerna-teplota-zari.png)
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
![](/images/c/ce/ZOD-cv12-prumerna-teplota-zari-pocet.png)
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
![](/images/e/ea/ZOD-cv12-prumerna-teplota-zari-filtr.png)
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
Souhrné informace o časové topologii vystiskne příkaz t.topology:
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 ...
Pro podrobné informace je nutné přidat přepínač -m:
t.topology -m modis where="start_time >= '2002-07-15' and start_time < '2002-07-17'"
+-------------------- Raster Dataset ----------------------------------------+ | | +-------------------- Basic information -------------------------------------+ | Id: ........................ aqua_lst_night20020715@modis2002lst_yfsg | Name: ...................... aqua_lst_night20020715 | Mapset: .................... modis2002lst_yfsg | Creator: ................... martin | Temporal type: ............. absolute | Creation time: ............. 2014-04-29 18:34:29.485248 +-------------------- Absolute time -----------------------------------------+ | Start time:................. 2002-07-15 01:30:00 | End time:................... None +-------------------- Temporal Topology -------------------------------------+ | Next: ...................... terra_lst_day20020715@modis2002lst_yfsg +-------------------- Spatial extent ----------------------------------------+ | North:...................... 279073.975466 | South:...................... 113673.975466 | East:.. .................... 798143.311797 | West:....................... 595143.311797 | Top:........................ 0.0 | Bottom:..................... 0.0 +-------------------- Metadata information ----------------------------------+ | Datatype:................... DCELL | Number of columns:.......... 827 | Number of rows:............. 1015 | Number of cells:............ 839405 | North-South resolution:..... 200.0 | East-west resolution:....... 200.0 | Minimum value:.............. 12.85 | Maximum value:.............. 19.53 | Registered datasets ........ modis@modis2002lst_yfsg +----------------------------------------------------------------------------+ ... +-------------------- Raster Dataset ----------------------------------------+ | | +-------------------- Basic information -------------------------------------+ | Id: ........................ terra_lst_night20020716@modis2002lst_yfsg | Name: ...................... terra_lst_night20020716 | Mapset: .................... modis2002lst_yfsg | Creator: ................... martin | Temporal type: ............. absolute | Creation time: ............. 2014-04-29 18:34:29.603983 +-------------------- Absolute time -----------------------------------------+ | Start time:................. 2002-07-16 22:30:00 | End time:................... None +-------------------- Temporal Topology -------------------------------------+ | Previous: .................. aqua_lst_day20020716@modis2002lst_yfsg +-------------------- Spatial extent ----------------------------------------+ | North:...................... 279073.975466 | South:...................... 113673.975466 | East:.. .................... 798143.311797 | West:....................... 595143.311797 | Top:........................ 0.0 | Bottom:..................... 0.0 +-------------------- Metadata information ----------------------------------+ | Datatype:................... DCELL | Number of columns:.......... 827 | Number of rows:............. 1015 | Number of cells:............ 839405 | North-South resolution:..... 200.0 | East-west resolution:....... 200.0 | Minimum value:.............. 18.95 | Maximum value:.............. 24.33 | Registered datasets ........ modis@modis2002lst_yfsg +----------------------------------------------------------------------------+
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|...
Agregace
Určení statististiky teplot pro jednotlivé měsíce:
t.rast.aggregate input=modis out=modis_m basename=ag gra="1 months"
Vytvoří se dvanáct rastrových map v měsíční periodě:
t.rast.list modis_m order=start_time
ag_0 modis2002lst_1 2002-01-01 00:00:00 2002-02-01 00:00:00 ag_1 modis2002lst_1 2002-02-01 00:00:00 2002-03-01 00:00:00 ag_2 modis2002lst_1 2002-03-01 00:00:00 2002-04-01 00:00:00 ag_3 modis2002lst_1 2002-04-01 00:00:00 2002-05-01 00:00:00 ag_4 modis2002lst_1 2002-05-01 00:00:00 2002-06-01 00:00:00 ag_5 modis2002lst_1 2002-06-01 00:00:00 2002-07-01 00:00:00 ag_6 modis2002lst_1 2002-07-01 00:00:00 2002-08-01 00:00:00 ag_7 modis2002lst_1 2002-08-01 00:00:00 2002-09-01 00:00:00 ag_8 modis2002lst_1 2002-09-01 00:00:00 2002-10-01 00:00:00 ag_9 modis2002lst_1 2002-10-01 00:00:00 2002-11-01 00:00:00 ag_10 modis2002lst_1 2002-11-01 00:00:00 2002-12-01 00:00:00 ag_11 modis2002lst_1 2002-12-01 00:00:00 2003-01-01 00:00:00
Statistika pro všechny měsíce:
t.rast.univar modis_m
ag_0@modis2002lst_1|2002-01-01 00:00:00|2002-02-01 00:00:00|10.1761589696136|6.11111111111114|14.3009090909091|... ag_11@modis2002lst_1|2002-12-01 00:00:00|2003-01-01 00:00:00|7.50005058325719|4.3941666666667|10.3609523809524|...
Statistika pouze pro červenec a srpen:
t.rast.univar modis_m where="start_time > '2002-07-01' and start_time < '2002-09-01'"
ag_6@modis2002lst_1|2002-07-01 00:00:00|2002-08-01 00:00:00|26.3466550685411|21.7652941176471|30.594|... ag_7@modis2002lst_1|2002-08-01 00:00:00|2002-09-01 00:00:00|26.7082631007969|22.2151851851852|31.4769565217392|...
Výběr dat z časoprostorového datasetu
Vytvořit na základě výběru nový časoprostorový dataset umožňuje příkaz t.rast.extract.
t.rast.extract input=modis where="start_time > '2002-03-01' and start_time < '2002-06-01'" output=modis_spring
t.rast.extract input=modis where="start_time > '2002-06-01' and start_time < '2002-09-01'" output=modis_summer
t.rast.extract input=modis where="start_time > '2002-09-01' and start_time < '2002-12-01'" output=modis_autumn
t.rast.extract input=modis where="start_time > '2002-12-01' or start_time < '2002-03-01'" output=modis_winter
Vizualizace časové řady pomocí nástroje g.gui.timeline.
g.gui.timeline inputs=modis_spring,modis_summer,modis_autumn,modis_winter
![](/images/7/7c/G-gui-timeline-modis-4.png)
V následujících příkazech budeme sledovat trend změny teploty v jednotlivých ročních obdobích. K tomu použijeme modul t.rast.series.
t.rast.series input=modis_spring output=modis_spring_avg method=average
t.rast.series input=modis_summer output=modis_summer_avg method=average
t.rast.series input=modis_autumn output=modis_autumn_avg method=average
t.rast.series input=modis_winter output=modis_winter_avg method=average
Vzniknou čtyři rastrové mapy zobrazující průměrné roční teploty v ročních obdobích. Průměrnou teplotu zjistíme pomocí modulu r.univar, příklad pro jaro:
r.univar modis_spring_avg
...
mean: 21.6125
...
![](/images/c/cf/Modis-spring-avg.png)
Vizualizace časoprostorových dat
Vizualizovat data časoprostorových datasetů umožňuje wxGUI animační nástroj.
g.gui.animation strds=modis
![](/images/e/ee/G-gui-animation.png)
Mezi další užitečné nástroje patří wxGUI Map Swipe
t.rast.list modis_m where="start_time < '2002-03-01'"
ag_0 modis2002lst_yfsg 2002-01-01 00:00:00 2002-02-01 00:00:00 ag_1 modis2002lst_yfsg 2002-02-01 00:00:00 2002-03-01 00:00:00
g.gui.mapswipe first=ag_0 second=ag_1
![](/images/4/43/G-gui-mapswipe-modis.png)
![](/images/7/70/G-gui-timeline-modis.png)