Landsat¶
V této části si ukážeme především proces importu reálných dat do systému GRASS, vytvoření časových značek a na jejich základě vytvoření časoprostorového datasetu.
Jako vstupní data máme sadu snímků Landsat ve formátu TIFF. Ke každému snímku je k dispozici maska, která definuje validního hodnoty. Maska obsahuje hodnoty 1 a 2, přičemž hodnota 2 představuje oblačnost. Např.
gdalinfo lndcal.LT51920262011306KIS00_20111102.tif -mm
gdalinfo fmask.LT51920262011306KIS00_20111102.tif -noct -mm
Viz Landsat Surface Reflectance Quality Assessment.
Tip
Existují užitečné utility, které usnadňují automatizované stažení dat Landsat. Patří mezi ně např. landsat-util, gsutil nebo Landsat-Download.
Import vstupních dat¶
V prvním kroku na základě vstupních dat vytvoříme novou lokaci:
grass70 -c lndcal.LT51920262011306KIS00_20111102.tif /opt/grassdata/landsat
Před importem můžeme zkontrolovat souřadnicový systém lokace:
g.proj -p
Pro dávkový import a vytvoření časoprostorového datasetu si napíšeme jednoduchý skript v jazyku Python, který provede následující kroky:
naimportuje data včetně masek 18
aplikuje masky nad originalními daty 34
nastaví časovou značku (datum je uvedeno jako součást názvu souboru) 59
maskované rastry zaregistruje do výstupního časoprostorového datasetu 73
1#!/usr/bin/env python
2
3#%module
4#% description: Creates temporal dataset from Landsat data.
5#%end
6#%option G_OPT_M_DIR
7#% description: Name of directory with input Landsat files
8#%end
9#%option G_OPT_STRDS_OUTPUT
10#%end
11
12import os
13import sys
14import datetime
15import grass.script as grass
16from grass.exceptions import CalledModuleError
17
18def import_data(directory):
19 files = []
20 for f in os.listdir(directory):
21 if f.endswith(".tif"):
22 files.append(f)
23
24 count = len(files)
25 i = 0
26 for f in files:
27 i += 1
28 grass.message("Importing <{0}> ({1}/{2})...".format(f, i, count))
29 grass.percent(i, count, 5)
30 grass.run_command('r.in.gdal', input=os.path.join(directory, f),
31 output=os.path.splitext(f)[0], quiet=True, overwrite=True)
32 grass.percent(0, 0, 0)
33
34def mask_data():
35 bands1 = grass.list_grouped('raster', pattern='lndcal*.1$')[mapset]
36 count = len(bands1)
37 i = 0
38 for b1 in bands1:
39 i += 1
40 oname = b1.split('.')[1]
41 grass.message("Processing <{0}> ({1}/{2})...".format(oname, i, count))
42 grass.percent(i, count, 5)
43 grass.run_command('g.region', raster=b1)
44 mask = grass.find_file('fmask.' + oname, element='raster')['fullname']
45 if mask:
46 grass.run_command('r.mask', raster='fmask.' + oname, maskcats=1, overwrite=True, quiet=True)
47 else:
48 grass.warning("Mask missing for <{0}>".format(oname))
49 bands = grass.list_grouped('raster', pattern='lndcal.{0}.*'.format(oname))[mapset]
50 for b in bands:
51 n = b.split('.')[-1]
52 grass.mapcalc('{0}.{1}={2}'.format(oname, n, b), quiet=True, overwrite=True)
53 if mask:
54 grass.run_command('r.mask', flags='r', quiet=True)
55
56 grass.run_command('g.remove', type='raster', pattern='fmask*,lndcal*', flags='f', quiet=True)
57 grass.percent(0, 0, 0)
58
59def timestamp_data():
60 maps = grass.list_grouped('raster')[mapset]
61 grass.message("Timestamping maps...")
62 count = len(maps)
63 i = 0
64 for name in maps:
65 i += 1
66 grass.percent(i, count, 5)
67 date = name.split('_')[-1].split('.')[0]
68 grass_date = datetime.datetime.strptime(date, '%Y%m%d').strftime('%d %b %Y')
69 grass.run_command('r.timestamp', map=name, date=grass_date)
70
71 grass.percent(0, 0, 0)
72
73def create_strds(output):
74 grass.run_command('t.create', output=output, title=output, description=output)
75 map_file = grass.tempfile()
76 grass.run_command('g.list', type='raster', mapset='.', separator='newline', output=map_file, overwrite=True)
77 grass.run_command('t.register', input=output, file=map_file, separator='newline', quiet=True, overwrite=True)
78
79def check_strds(name):
80 strds = grass.read_command('t.list', quiet=True).splitlines()
81 if name + '@' + mapset not in strds:
82 return
83
84 if os.environ.get('GRASS_OVERWRITE', '0') == '1':
85 grass.warning("Space time raster dataset <{}> is already exists "
86 "and will be overwritten.".format(name))
87 grass.run_command('t.remove', inputs=name, quiet=True)
88 else:
89 grass.fatal("Space time raster dataset <{}> is already in the database. "
90 "Use the overwrite flag.".format(name))
91
92def main():
93 check_strds(opt['output'])
94 import_data(opt['input'])
95 mask_data()
96 timestamp_data()
97 create_strds(opt['output'])
98
99if __name__ == "__main__":
100 opt, flgs = grass.parser()
101 mapset=grass.gisenv()['MAPSET']
102 sys.exit(main())
Skript je ke stažení zde. Příklad spuštění:
grass_landsat.py input=~/geodata/sub101 out=landsat
Ve výsledném datasetu máme zaregistrována data od roku 1984 do 2014:
t.info landsat
...
| Start time:................. 1984-06-25 00:00:00
| End time:................... 2014-12-29 00:00:00
...
Tip
Vylepšená verze skriptu s podporou pro multiprocessing je ke stažení zde.
Detektce vodní ploch v časové řadě¶
t.rast.extract input=landsat where="name like '%_B2'" output=landsat_b2
t.rast.extract input=landsat where="name like '%_B5'" output=landsat_b5
t.rast.mapcalc inputs=landsat_b2,landsat_b5
expression="float(landsat_b2 - landsat_b5) / (landsat_b2 + landsat_b5) > 0"
output=landsat_voda basename=voda