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
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
#!/usr/bin/env python

#%module
#% description: Creates temporal dataset from Landsat data.
#%end
#%option G_OPT_M_DIR
#% description: Name of directory with input Landsat files
#%end
#%option G_OPT_STRDS_OUTPUT
#%end

import os
import sys
import datetime
import grass.script as grass
from grass.exceptions import CalledModuleError

def import_data(directory):
    files = []
    for f in os.listdir(directory):
        if f.endswith(".tif"):
            files.append(f)

    count = len(files)
    i = 0
    for f in files:
        i += 1
        grass.message("Importing <{0}> ({1}/{2})...".format(f, i, count))
        grass.percent(i, count, 5)
        grass.run_command('r.in.gdal', input=os.path.join(directory, f),
                          output=os.path.splitext(f)[0], quiet=True, overwrite=True)
    grass.percent(0, 0, 0)
    
def mask_data():
    bands1 = grass.list_grouped('raster', pattern='lndcal*.1$')[mapset]
    count = len(bands1)
    i = 0
    for b1 in bands1:
        i += 1
        oname = b1.split('.')[1]
        grass.message("Processing <{0}> ({1}/{2})...".format(oname, i, count))
        grass.percent(i, count, 5)
        grass.run_command('g.region', raster=b1)
        mask = grass.find_file('fmask.' + oname, element='raster')['fullname']
        if mask:
            grass.run_command('r.mask', raster='fmask.' + oname, maskcats=1, overwrite=True, quiet=True)
        else:
            grass.warning("Mask missing for <{0}>".format(oname))
        bands = grass.list_grouped('raster', pattern='lndcal.{0}.*'.format(oname))[mapset]
        for b in bands:
            n = b.split('.')[-1]
            grass.mapcalc('{0}.{1}={2}'.format(oname, n, b), quiet=True, overwrite=True)
        if mask:
            grass.run_command('r.mask', flags='r', quiet=True)
        
    grass.run_command('g.remove', type='raster', pattern='fmask*,lndcal*', flags='f', quiet=True)
    grass.percent(0, 0, 0)
    
def timestamp_data():
    maps = grass.list_grouped('raster')[mapset]
    grass.message("Timestamping maps...")
    count = len(maps)
    i = 0
    for name in maps:
        i += 1
        grass.percent(i, count, 5)
        date = name.split('_')[-1].split('.')[0]
        grass_date = datetime.datetime.strptime(date, '%Y%m%d').strftime('%d %b %Y')
        grass.run_command('r.timestamp', map=name, date=grass_date)

    grass.percent(0, 0, 0)

def create_strds(output):
    grass.run_command('t.create', output=output, title=output, description=output)
    map_file = grass.tempfile()
    grass.run_command('g.list', type='raster', mapset='.', separator='newline', output=map_file, overwrite=True)
    grass.run_command('t.register', input=output, file=map_file, separator='newline', quiet=True, overwrite=True)

def check_strds(name):
    strds = grass.read_command('t.list', quiet=True).splitlines()
    if name + '@' + mapset not in strds:
        return
    
    if os.environ.get('GRASS_OVERWRITE', '0') == '1':
        grass.warning("Space time raster dataset <{}> is already exists "
                      "and will be overwritten.".format(name))
        grass.run_command('t.remove', inputs=name, quiet=True)
    else:
        grass.fatal("Space time raster dataset <{}> is already in the database. "
                    "Use the overwrite flag.".format(name))
    
def main():
    check_strds(opt['output'])
    import_data(opt['input'])
    mask_data()
    timestamp_data()
    create_strds(opt['output'])
    
if __name__ == "__main__":
    opt, flgs = grass.parser()
    mapset=grass.gisenv()['MAPSET']
    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