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:

  1. naimportuje data včetně masek 18

  2. aplikuje masky nad originalními daty 34

  3. nastaví časovou značku (datum je uvedeno jako součást názvu souboru) 59

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