Difference between revisions of "LANDSAT"

From GRASS-Wiki
Jump to: navigation, search
(Modules)
(Notes)
Line 38: Line 38:
  
 
The {{AddonCmd|i.landsat.toar}} and {{AddonCmd|i.landsat.acca}} modules want the maps to be named such as "scenename.1, .2, .3", etc. for the different bands. GloVis names LANDSAT-7 like _B10, _B20, _B30, etc.
 
The {{AddonCmd|i.landsat.toar}} and {{AddonCmd|i.landsat.acca}} modules want the maps to be named such as "scenename.1, .2, .3", etc. for the different bands. GloVis names LANDSAT-7 like _B10, _B20, _B30, etc.
 +
 +
=== Import data ===
 +
 +
An example of Python script bellow imports raster data into GRASS. For each image creates its own mapset and imports bands as <tt>B<id></tt>, e.g. B10, B20, ..., B80. This script also sets up timestamp based on MTL file.
 +
 +
<source lang=python>
 +
#!/usr/bin/python
 +
 +
import os
 +
import sys
 +
import glob
 +
import grass.script as grass
 +
 +
def get_timestamp(mapset):
 +
    try:
 +
        metafile = glob.glob(mapset + '/*MTL.txt')[0]
 +
    except IndexError:
 +
        return
 +
 +
    result = dict()
 +
    try:
 +
        fd = open(metafile)
 +
        for line in fd.readlines():
 +
            line = line.rstrip('\n')
 +
            if len(line) == 0:
 +
                continue
 +
            if 'ACQUISITION_DATE' in line:
 +
                result['date'] = line.strip().split('=')[1].strip()
 +
    finally:
 +
        fd.close()
 +
 +
    return result
 +
 +
def import_tifs(mapset):
 +
    meta = get_timestamp(mapset)
 +
    for file in os.listdir(mapset):
 +
        if os.path.splitext(file)[-1] != '.TIF':
 +
            continue
 +
        ffile = os.path.join(mapset, file)
 +
        name = os.path.splitext(file)[0].split('_')[-1]
 +
        kanal = int(name[-2])
 +
        grass.message('Importuji %s -> %s@%s...' % (file, name, mapset))
 +
        grass.run_command('g.mapset',
 +
                          flags = 'c',
 +
                          mapset = mapset,
 +
                          quiet = True,
 +
                          stderr = open(os.devnull, 'w'))
 +
        grass.run_command('r.in.gdal',
 +
                          input = ffile,
 +
                          output = name,
 +
                          quiet = True,
 +
                          overwrite = True,
 +
                          title = 'kanal %d' % kanal)
 +
        if meta:
 +
            year, month, day = meta['date'].split('-')
 +
            if month == '01':
 +
                month = 'jan'
 +
            elif month == '02':
 +
                month = 'feb'
 +
            elif month == '03':
 +
                month = 'mar'
 +
            elif month == '04':
 +
                month = 'apr'
 +
            elif month == '05':
 +
                month = 'may'
 +
            elif month == '06':
 +
                month = 'jun'
 +
            elif month == '07':
 +
                month = 'jul'
 +
            elif month == '08':
 +
                month = 'aug'
 +
            elif month == '09':
 +
                month = 'sep'
 +
            elif month == '10':
 +
                month = 'oct'
 +
            elif month == '11':
 +
                month = 'nov'
 +
            elif month == '12':
 +
                month = 'dec'
 +
 +
            grass.run_command('r.timestamp',
 +
                              map = name,
 +
                              date = ' '.join((day, month, year)))
 +
 +
def main():
 +
    if len(sys.argv) == 1:
 +
        for directory in filter(os.path.isdir, os.listdir(os.getcwd())):
 +
            import_tifs(directory)
 +
    else:
 +
        import_tifs(sys.argv[1])
 +
 +
if __name__ == "__main__":
 +
    main()
 +
</source>
 +
 +
Example of usage:
 +
 +
* <tt>./import.py</tt> walk through current directory and import all found satellite images.
 +
* <tt>./import.py LM41890261983200FFF03</tt> imports images only from given directory.
  
 
=== Minimal disk space copies ===
 
=== Minimal disk space copies ===

Revision as of 06:55, 24 October 2010

LANDSAT Data Availability

Modules

  • d.rgb - display 3-band data
  • i.landsat.rgb - auto-enhance colors
  • i.atcorr - correct top of atmosphere to surface reflectance
  • i.oif - calculate the 3 bands showing the greatest difference (for use as R,G,B bands)
  • r.composite - flatten 3-bands of data into a single image (lossy)

LANDSAT Pre-Processing

Some notes from Yann Chemin:

  1. Open GRASS GIS, select create location from georeferenced file
  2. Use r.in.gdal to import into GRASS GIS.
    1. Most landsat scenes are delivered in "north-is-up" orientation, hence import is straightforward
    2. If you get "ERROR: Input map is rotated - cannot import.", the image must be first externally rotated, applying the rotation info stored in the metadata field of the raster image file. For example, the gdalwarp software can be used to transform the map to North-up (note, there are several gdalwarp parameters to select the resampling algorithm):
 gdalwarp rotated.tif northup.tif
  1. In grass/add-ons look for any of these to correct from radiance to reflectance at top of atmosphere
    1. i.landsat.toar
    2. i.dn2full.l5
    3. i.dn2full.l7
  2. Use i.atcorr to correct top of atmosphere to surface reflectance

Notes

The i.landsat.toar and i.landsat.acca modules want the maps to be named such as "scenename.1, .2, .3", etc. for the different bands. GloVis names LANDSAT-7 like _B10, _B20, _B30, etc.

Import data

An example of Python script bellow imports raster data into GRASS. For each image creates its own mapset and imports bands as B<id>, e.g. B10, B20, ..., B80. This script also sets up timestamp based on MTL file.

#!/usr/bin/python
 
import os
import sys
import glob
import grass.script as grass
 
def get_timestamp(mapset):
    try:
        metafile = glob.glob(mapset + '/*MTL.txt')[0]
    except IndexError:
        return
 
    result = dict()
    try:
        fd = open(metafile)
        for line in fd.readlines():
            line = line.rstrip('\n')
            if len(line) == 0:
                continue
            if 'ACQUISITION_DATE' in line:
                result['date'] = line.strip().split('=')[1].strip()
    finally:
        fd.close()
 
    return result
 
def import_tifs(mapset):
    meta = get_timestamp(mapset)
    for file in os.listdir(mapset):
        if os.path.splitext(file)[-1] != '.TIF':
            continue
        ffile = os.path.join(mapset, file)
        name = os.path.splitext(file)[0].split('_')[-1]
        kanal = int(name[-2])
        grass.message('Importuji %s -> %s@%s...' % (file, name, mapset))
        grass.run_command('g.mapset',
                          flags = 'c',
                          mapset = mapset,
                          quiet = True,
                          stderr = open(os.devnull, 'w'))
        grass.run_command('r.in.gdal',
                          input = ffile,
                          output = name,
                          quiet = True,
                          overwrite = True,
                          title = 'kanal %d' % kanal)
        if meta:
            year, month, day = meta['date'].split('-')
            if month == '01':
                month = 'jan'
            elif month == '02':
                month = 'feb'
            elif month == '03':
                month = 'mar'
            elif month == '04':
                month = 'apr'
            elif month == '05':
                month = 'may'
            elif month == '06':
                month = 'jun'
            elif month == '07':
                month = 'jul'
            elif month == '08':
                month = 'aug'
            elif month == '09':
                month = 'sep'
            elif month == '10':
                month = 'oct'
            elif month == '11':
                month = 'nov'
            elif month == '12':
                month = 'dec'
 
            grass.run_command('r.timestamp',
                              map = name,
                              date = ' '.join((day, month, year)))
 
def main():
    if len(sys.argv) == 1:
        for directory in filter(os.path.isdir, os.listdir(os.getcwd())):
            import_tifs(directory)
    else:
        import_tifs(sys.argv[1])
 
if __name__ == "__main__":
    main()

Example of usage:

  • ./import.py walk through current directory and import all found satellite images.
  • ./import.py LM41890261983200FFF03 imports images only from given directory.

Minimal disk space copies

Here's a little trick with r.reclass to rename it without touching the data or wasting disk space:

BASE=L71074092_09220040924

for BAND in 10 20 30 40 50 61 70 80; do
  BAND1st=`echo $BAND | sed -e 's/0$//'`
  r.reclass in="${BASE}_B$BAND" out=$BASE.$BAND1st << EOF
    * = *
EOF
done


Reset color tables

for map in `g.mlist pat="$BASE.[0-8]*"` ; do
  r.colors $map color=grey255
done


Create an MASK to only show data where coverage exists for all bands

g.region rast=$BASE.1
r.series in=`g.mlist pat="$BASE.[0-8]*" sep=,` -n out=$BASE.thresh method=threshold thresh=1
r.mapcalc "$BASE.mask = if(isnull($BASE.thresh))"
g.remove $BASE.thresh



Calculate top-of-atmosphere reflectance and band-6 temperature:

i.landsat.toar sensor=7 band_pre=$BASE metfile=${BASE}_MTL.txt

Identify clouds in the image:

i.landsat.acca -s -f band_prefix=$BASE.toar out=$BASE.acca

Mask out the clouds:

r.mapcalc "MASK = if(isnull($BASE.acca))"

Sample data

The North Carolina 2008 sample dataset comes with 3 different Landsat scenes:

(Wake County -- path: 16 row: 35 for various dates)

Landsat-5: Oct 14, 1987

Glovis download
LT50160351987287XXX08
http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LT50160351987287XXX08&dataset_name=LANDSAT_TM


Provided metadata:

IMAGE_ID=P016R35_5T871014
PATH=16
ROW=35
DATE=10/14/87
PLATFORM=LANDSAT5

Landsat-7: Mar 31, 2000

Glovis download
LE70160352000091EDC00
http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352000091EDC00&dataset_name=LANDSAT_ETM

Values in the metadata below indicate that the version provided with the NC 2008 dataset's production date was after July 1, 2000, and that the channel gains were HHHLHLHHL.

Convert DNs to radiance/temperatures:

GRASS> i.landsat.toar -v band=lsat7_2000 sensor=7 date=2000-03-31 \
   product_date=2000-07-02 solar_elevation=51.5246529 gain=HHHLHLHHL

Provided metadata:

SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
ACQUISITION_DATE=2000-03-31
WRS_PATH=16
CPF_FILE_NAME=L7CPF20000101_20000331_12
SUN_AZIMUTH=139.6033279
SUN_ELEVATION=51.5246529
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0

Landsat-7: May 24, 2002

Glovis download
LE70160352002144EDC00
http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352002144EDC00&dataset_name=LANDSAT_ETM


Provided metadata: (`p016r035_7x20020524.met`)

SPACECRAFT_ID=Landsat7
SENSOR_ID=ETM+
ACQUISITION_DATE=2002-05-24
WRS_PATH=016
WRS_ROW=035
SUN_AZIMUTH=120.8810347
SUN_ELEVATION=64.7730999
QA_PERCENT_MISSING_DATA=0
CLOUD_COVER=0
CPF_FILE_NAME=L7CPF20020401_20020630_03
LMAX_BAND1=191.600
LMIN_BAND1=-6.200
LMAX_BAND2=196.500
LMIN_BAND2=-6.400
LMAX_BAND3=152.900
LMIN_BAND3=-5.000
LMAX_BAND4=241.100
LMIN_BAND4=-5.100
LMAX_BAND5=31.060
LMIN_BAND5=-1.000
LMAX_BAND61=17.040
LMIN_BAND61=0.000
LMAX_BAND62=12.650
LMIN_BAND62=3.200
LMAX_BAND7=10.800
LMIN_BAND7=-0.350
LMAX_BAND8=243.100
LMIN_BAND8=-4.700
QCALMAX_BAND1=255.0
QCALMIN_BAND1=1.0
QCALMAX_BAND2=255.0
QCALMIN_BAND2=1.0
QCALMAX_BAND3=255.0
QCALMIN_BAND3=1.0
QCALMAX_BAND4=255.0
QCALMIN_BAND4=1.0
QCALMAX_BAND5=255.0
QCALMIN_BAND5=1.0
QCALMAX_BAND61=255.0
QCALMIN_BAND61=1.0
QCALMAX_BAND62=255.0
QCALMIN_BAND62=1.0
QCALMAX_BAND7=255.0
QCALMIN_BAND7=1.0
QCALMAX_BAND8=255.0
QCALMIN_BAND8=1.0