LANDSAT: Difference between revisions

From GRASS-Wiki
Jump to navigation Jump to search
mNo edit summary
 
(100 intermediate revisions by 9 users not shown)
Line 1: Line 1:
== LANDSAT Data Availability==
== Data Availability==


* data download from [http://landsat.usgs.gov/ USGS LANDSAT archive] and [http://glcfapp.umiacs.umd.edu GLCF]
Landsat imagery can be obtained from [http://landsat.usgs.gov/ USGS' LANDSAT archive]:
* USGS GloVis Java app: http://glovis.usgs.gov
* see also the [[Global datasets]] wiki page


== Modules ==
* [http://glovis.usgs.gov USGS GloVis]
* [http://earthexplorer.usgs.gov/ EarthExplorer]
* [http://earthexplorer.usgs.gov/bulk/help Bulk Download Orders] using the [http://earthexplorer.usgs.gov/bulk/ Bulk Download Application]
* [http://landsatlook.usgs.gov/ LandsatLook Viewer]
* [https://landsatonaws.com/ Landsat on Amazon Web Services]


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


* {{AddonCmd|i.landsat.toar}} (addon) - convert DN to top of atmosphere radiance
* the [[Global datasets]] wiki page
* {{AddonCmd|i.landsat.acca}} (addon) - cloud identification
* [http://landsat.usgs.gov/band_designations_landsat_satellites.php Landsat satellites band designations].
* {{AddonCmd|i.landsat.dehaze}} (addon) - haze removal


== LANDSAT Pre-Processing ==
== Modules overview ==


Some notes from Yann Chemin:
GRASS GIS features a complete set of native modules and various add-ons for pre- and post-processing Landsat satellite imagery. The following lists offer an overview of related modules and add-ons.


# Open GRASS GIS, select create location from georeferenced file
=== Generic modules applicable to Landsat ===
# Use {{cmd|r.in.gdal}} to import into GRASS GIS.
 
## Most landsat scenes are delivered in "north-is-up" orientation, hence import is straightforward
* {{cmd|d.rgb}} → display 3-band data
## 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):
* {{cmd|r.composite}} → flatten 3-bands of data into a single image (lossy, maybe used in combination with {{cmd|i.colors.enhance}}
  gdalwarp rotated.tif northup.tif
* {{cmd|i.atcorr}} → correct top of atmosphere to surface reflectance - see also the [[Atmospheric correction]] wiki page
* {{cmd|i.topo.corr}} → used to topographically correct reflectance from imagery files, e.g. obtained with {{cmd|i.landsat.toar}}, using a sun illumination terrain model
* {{cmd|i.oif}} → calculate the 3 bands showing the greatest difference (for use as R,G,B bands)
* {{AddonCmd|i.histo.match}} (addon) → histogram matching of two or more raster maps (in grass 7) ''Note, the module works with integer values and does not accept the "." character as part of the raster map's name!''
 
=== Landsat specific modules ===
 
* {{cmd|i.colors.enhance}} → auto-enhance colors
* {{cmd|i.landsat.toar}} → convert DN to top of atmosphere radiance
* {{cmd|i.landsat.acca}} → cloud cover assessment
* {{cmd|i.tasscap}} → Performs Tasseled Cap (Kauth Thomas) transformation
 
=== Landsat specific GRASS AddOns ===
 
* {{AddonCmd|i.landsat.trim|version=7}} → trims border fringes for each band separately or with the MASK where coverage exists for all bands
* {{AddonCmd|i.landsat.dehaze|version=7}} → haze removal
* {{AddonCmd|i.landsat8.qc|version=7}} → Reclass Landsat8 QA band according to acceptable pixel quality as defined by the user
* {{AddonCmd|i.landsat8.swlst|version=7}} → Practical split-window algorithm estimating Land Surface Temperature from Landsat 8 OLI/TIRS imagery
* {{AddonCmd|i.wi|version=7}} → Calculate different water indices
 
== Pre-Processing ==
 
=== Overview ===
 
Typically, pre-processing Landsat imagery comprises the following steps:
 
# '''import''' in the database → {{cmd|r.in.gdal}}
# geometrically & orthometrically correct imagery
#* already done for L1T products, read more at [http://landsat.usgs.gov//Landsat_Processing_Details.php USGS' Landsat Information Products] webpage
# optionally, automatically '''cut-off border fringes''' → {{AddonCmd|i.landsat.trim|version=7}}
#* of course one can use the official WRS2 Path/Row vector tiles to manually trim border fringes
# optionally, '''denoise''' for obvious/intensive salt & pepper effects, stripes, etc.
#* for example by applying Principal Components Analysis as a denoising technique → {{cmd|i.pca}} ''<<< Re-order this step?''
# '''convert''' the '''Digital Numbers''' (DNs) to '''Top-of-Atmosphere Radiances/Reflectances''' (ToARs) → {{cmd|i.landsat.toar}}
# optionally, '''correct for atmospheric effects''' → {{cmd|i.atcorr}}
#* that is, accounting for distorting atmospheric effects and estimating actual reflectances as they would have been measured on the ground
#* also described as conversion to Top-of-Canopy Reflectances (ToCRs)
# '''assessing cloud cover''' → {{cmd|i.landsat.acca}}
#* optionally, detect and remove clouds shadows as well
# '''topographically normalise''' imagery → {{cmd|i.topo.corr}}
#* also known as topographic correction, that is, accounting for illumination differences due to the acquisition's geometry
# '''radiometrically normalise''' → one approach via {{AddonCmd|i.histo.match}} (in '''grass 7''', addon), also known as relative radiometric normalisation -- one approach is the ''histogram matching'' technique of two or more raster maps
#* '''ToDo/More techniques?'''
 
See also
* [http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/ Processing Landsat 8 data in GRASS GIS 7: Import and visualization]
 
=== Importing data ===
 
Importing Landsat spectral bands in GRASS GIS' data base can be done both from the Graphical User Interface or from the command line
 
 
# Open GRASS GIS, select ''Location Wizard'' in order to create a new location from georeferenced file
# Use the ''Import file tool'', or the {{cmd|r.in.gdal}} module, to import the GeoTIFF files into GRASS GIS.
 
See also [http://grasswiki.osgeo.org/wiki/LANDSAT#Automated_data_import Automated data import] below.


# In grass/add-ons look for any of these to correct from radiance to reflectance at top of atmosphere
==== Notes ====
## {{AddonCmd|i.landsat.toar}}
## {{AddonCmd|i.dn2full.l5}}
## {{AddonCmd|i.dn2full.l7}}
# Use {{cmd|i.atcorr}} to correct top of atmosphere to surface reflectance


== Notes ==
# Most Landsat scenes are delivered in ''north-is-up'' orientation, hence the import process is straightforward.
# If you get <source lang="bash" enclose="none">ERROR: Input map is rotated - cannot import.</source>, the image must be first rotated externally, applying the rotation info stored in the metadata field of the raster image file. For example, the <code>gdalwarp</code> software can be used to transform the map to North-up (note, there are several <code>gdalwarp</code> parameters to select the resampling algorithm):


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.
<ul>
<source lang="bash">
  gdalwarp rotated.tif northup.tif
</source>
</ul>


Here's a little trick with {{cmd|r.reclass}} to rename it without touching the data or wasting disk space:
==== Hint: Minimal disk space copies ====
 
Here's a little trick using {{cmd|r.reclass}} to rename maps (for example, from ''L71074092_09220040924_B'''10''''' to ''L71074092_09220040924_B'''1''''') without touching the data or wasting disk space:


<source lang="bash">
<source lang="bash">
Line 50: Line 105:
</source>
</source>


Calculate top-of-atmosphere reflectance and band-6 temperature:
i.landsat.toar -7 band_pre=$BASE metfile=${BASE}_MTL.txt


Identify clouds in the image:
 
  i.landsat.acca -2 -f band_prefix=$BASE.toar out=$BASE.acca
==== Automated data import ====
 
The following ''Python'' script imports Landsat imagery into GRASS' data base. Specifically, the script
 
* creates an independent Mapset for each Landsat scene
* imports and renames bands of a scene as <tt>B<id></tt>, e.g. B10, B20, ..., B80.
* additionaly sets up the timestamp based on MTL metadata file
 
 
'''Note,''' the (newest) official naming pattern for Landsat scenes -- explained in [https://lta.cr.usgs.gov/landsat_dictionary.html USGS' Landsat Data Dictionary] as the [https://lta.cr.usgs.gov/landsat_dictionary.html#entity_id Landsat Scene Identifier] -- and all individual bands that compose a scene -- bands have a suffix which is like _B10, _B20, _B30, etc. -- differ from what some Landsat specific modules expect. For example, the modules {{cmd|i.landsat.toar}} and {{cmd|i.landsat.acca}} expect the bands to follow a naming pattern such as "scenename.1, .2, .3", etc.
 
 
To use the script save it as <source lang="bash" enclose="none">import_landsat.py</source> file and make sure it is granted the execution permission.
 
<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 any(x in line for x in ('DATE_ACQUIRED', 'ACQUISITION_DATE')):
                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)
        if ('VCID') in ffile:
            name = "".join((os.path.splitext(file)[0].split('_'))[1::2])
        else:
            name = os.path.splitext(file)[0].split('_')[-1]
        if len(name) == 3 and name[-1] == '0':
            band = int(name[1:2])
        elif len(name) == 3 and name[-1] != '0':
            band = int(name[1:3])
        else:
            band = int(name[-1:])
        grass.message('Importing %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 = 'band %d' % band)
        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'''
 
After having collected the Landsat scenes of interest in one directory, the script can be used as follows:
 
* <source lang="python" enclose="none">./import_landsat.py</source> → the script will walk through a ''pool'' directory that contains unique Landsat scene directories (e.g. three directories named after the official naming pattern: LT51800342011158MOR00 LT51810352009079MTI00  LT51820352009326MTI00) and import all bands of each individual scene in their own Mapset
* <source lang="python" enclose="none">./import_landsat.py LM41890261983200FFF03</source> → the scrip will import only bands of the specified Landsat scene directory
 
== Post-Processing ==
 
=== Natural color composites ===
 
Creating natural (also known as true-) color composites, can be done by
 
a) using {{cmd|i.colors.enhance}} to automatically balance the colors for the Red, Green and Blue bands
 
or
 
b) equalizing the colors on each of the Red, Green and Blue bands
 
<ul>
<source lang="bash">
r.colors -e map=band1 color=grey
</source>
</ul>
 
'''Composites''' can then be produced with the {{cmd|r.composite}} module.
 
'''Resetting the color tables''' of all bands back to normal greyscale can be done with a for loop:
 
<source lang="bash">
BASE=L71074092_09220040924
 
for map in `g.mlist pat="$BASE.[0-8]*"` ; do
  r.colors $map color=grey255
done
</source>
 
See also:
* [http://courses.neteler.org/processing-landsat-8-data-in-grass-gis-7-rgb-composites-and-pan-sharpening/ Processing Landsat 8 data in GRASS GIS 7: RGB composites and pan sharpening]
 
=== Create a MASK to only show data where coverage exists for all bands ===
 
<source lang="bash">
BASE=L71074092_09220040924
 
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
</source>
 
=== Calculate Top-of-Atmosphere Reflectance and band-6 Temperature ===
 
Calculate ''Top-of-Atmosphere'' reflectance and band-6 temperature using the {{cmd|i.landsat.toar}} module. For details, refer to USGS' [http://landsat.usgs.gov/Landsat_Metadata_Changes.php Landsat Filename and Metadata Changes] dedicated page. For Landsat 8, see also [http://landsat.usgs.gov/about_LU_Vol_7_Issue_4.php#3a here].
<source lang="bash">
i.landsat.toar input_prefix=$BASE output_prefix=${BASE}_toar metfile=${BASE}_MTL.txt
</source>
 
Note, the resulting temperature map is in Kelvin:
<source lang="bash">
# convert to degree Celsius
r.mapcalc "$BASE.temp_celsius = ${BASE}_toar.6 - 273.15"
r.info -r $BASE.temp_celsius
</source>
 
=== Haze removal ===
 
Simple haze removal can be done with {{AddonCmd|i.landsat.dehaze|version=7}}. This addons applies a bandwise haze correction using tasscap4 (haze) and linear regression.
 
During the 2000s, prior acquired Landsat data were reprocessed to LPGS (Level 1 Product Generation System). Seems that with this level of processing, haze effect in raw data was removed and now is sufficient to apply {{Cmd|i.colors.enhance}} to obtain a good image with high contrast.
 
=== Atmospheric correction ===
 
See [[Atmospheric correction]]
 
=== Cloud identification ===
 
Identify clouds in the image with {{cmd|i.landsat.acca}}:
<source lang="bash">i.landsat.acca -f input_prefix=226_62_toar. output=226_62.acca_cloudmask</source>


Mask out the clouds:
Mask out the clouds:
r.mapcalc "MASK = if( isnull($BASE.acca), 1, null() )"
<source lang="bash">r.mapcalc "MASK = if(isnull($BASE.acca_cloudmask))"</source>


== Download sample data ==


== Sample data ==
=== Preprocessed Landsat-7 data for North Carolina ===


The [[Sample_datasets|North Carolina 2008 sample dataset]] comes with 3 different Landsat scenes:
The [[Sample_datasets|North Carolina 2008 sample dataset]] comes with 3 different Landsat scenes:
: (''Wake County -- path: 16  row: 35 for various dates'')
The above import efforts are not needed since the data are already in a GRASS location.


=== Landsat-5: Oct 14, 1987 ===
=== Landsat-5: Oct 14, 1987 ===


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


Provided metadata:
Provided metadata:
Line 77: Line 328:
=== Landsat-7: Mar 31, 2000 ===
=== Landsat-7: Mar 31, 2000 ===


;Glovis download:
;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 production date was after July 1, 2000, and that the channel gains were <tt>HHHLHLHHL</tt>.
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 <tt>HHHLHLHHL</tt>.


Convert DNs to radiance/temperatures:
Convert DNs to radiance/temperatures:<source lang="bash">
  GRASS> i.landsat.toar -v band=lsat7_2000 sensor=7 date=2000-03-31 \
  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
     product_date=2000-07-02 solar_elevation=51.5246529 gain=HHHLHLHHL</source>


Provided metadata:
Provided metadata:
Line 130: Line 382:
  QCALMIN_BAND8=1.0
  QCALMIN_BAND8=1.0


=== Landsat-7: Apr 24, 2002 ===
=== Landsat-7: May 24, 2002 ===


;Glovis download:
;Glovis download: LE70160352002144EDC00
: http://edcsns17.cr.usgs.gov/cgi-bin/EarthExplorer/run-phtml/results/download.phtml?node=GV&ordered=LE70160352002144EDC00&dataset_name=LANDSAT_ETM
 
<!-- ???
or ELP016R035_7T20020524
: http://edcsns17.cr.usgs.gov/helpdocs/download.cgi?file=lsatortho/etm/16/35/elp016r035_7t20020524.tar.gz
--->


Provided metadata: (`p016r035_7x20020524.met`)
Provided metadata: (`p016r035_7x20020524.met`)
Line 182: Line 440:
  QCALMIN_BAND8=1.0
  QCALMIN_BAND8=1.0


== LANDSAT Image classification ==
See [[Image classification]]
== Time series analysis ==
See [[Time series]]


== Sources ==


* [http://landsat.usgs.gov/tools_access_all_faqs.php Frequently Asked Questions about the Landsat Missions]
* [https://landsat8portal.eo.esa.int/portal/ European LANDSAT 8 data in near-realtime]
* [http://courses.neteler.org/processing-landsat8-data-in-grass-gis-7/ Processing Landsat 8 data in GRASS GIS 7: Import and visualization]
* [http://courses.neteler.org/processing-landsat-8-data-in-grass-gis-7-rgb-composites-and-pan-sharpening/ Processing Landsat 8 data in GRASS GIS 7: RGB composites and pan sharpening]


[[Category:Documentation]]
[[Category: Documentation]]
[[Category:FAQ]]
[[Category: FAQ]]
[[Category:Image processing]]
[[Category: Landsat]]
[[Category: Image processing]]
[[Category: Import]]
[[Category: Digital Numbers]]
[[Category: Radiance]]
[[Category: Reflectance]]
[[Category: Geodata]]

Latest revision as of 10:31, 4 December 2018

Data Availability

Landsat imagery can be obtained from USGS' LANDSAT archive:

See also

Modules overview

GRASS GIS features a complete set of native modules and various add-ons for pre- and post-processing Landsat satellite imagery. The following lists offer an overview of related modules and add-ons.

Generic modules applicable to Landsat

  • d.rgb → display 3-band data
  • r.composite → flatten 3-bands of data into a single image (lossy, maybe used in combination with i.colors.enhance
  • i.atcorr → correct top of atmosphere to surface reflectance - see also the Atmospheric correction wiki page
  • i.topo.corr → used to topographically correct reflectance from imagery files, e.g. obtained with i.landsat.toar, using a sun illumination terrain model
  • i.oif → calculate the 3 bands showing the greatest difference (for use as R,G,B bands)
  • i.histo.match (addon) → histogram matching of two or more raster maps (in grass 7) Note, the module works with integer values and does not accept the "." character as part of the raster map's name!

Landsat specific modules

Landsat specific GRASS AddOns

  • i.landsat.trim → trims border fringes for each band separately or with the MASK where coverage exists for all bands
  • i.landsat.dehaze → haze removal
  • i.landsat8.qc → Reclass Landsat8 QA band according to acceptable pixel quality as defined by the user
  • i.landsat8.swlst → Practical split-window algorithm estimating Land Surface Temperature from Landsat 8 OLI/TIRS imagery
  • i.wi → Calculate different water indices

Pre-Processing

Overview

Typically, pre-processing Landsat imagery comprises the following steps:

  1. import in the database → r.in.gdal
  2. geometrically & orthometrically correct imagery
  3. optionally, automatically cut-off border fringesi.landsat.trim
    • of course one can use the official WRS2 Path/Row vector tiles to manually trim border fringes
  4. optionally, denoise for obvious/intensive salt & pepper effects, stripes, etc.
    • for example by applying Principal Components Analysis as a denoising technique → i.pca <<< Re-order this step?
  5. convert the Digital Numbers (DNs) to Top-of-Atmosphere Radiances/Reflectances (ToARs) → i.landsat.toar
  6. optionally, correct for atmospheric effectsi.atcorr
    • that is, accounting for distorting atmospheric effects and estimating actual reflectances as they would have been measured on the ground
    • also described as conversion to Top-of-Canopy Reflectances (ToCRs)
  7. assessing cloud coveri.landsat.acca
    • optionally, detect and remove clouds shadows as well
  8. topographically normalise imagery → i.topo.corr
    • also known as topographic correction, that is, accounting for illumination differences due to the acquisition's geometry
  9. radiometrically normalise → one approach via i.histo.match (in grass 7, addon), also known as relative radiometric normalisation -- one approach is the histogram matching technique of two or more raster maps
    • ToDo/More techniques?

See also

Importing data

Importing Landsat spectral bands in GRASS GIS' data base can be done both from the Graphical User Interface or from the command line


  1. Open GRASS GIS, select Location Wizard in order to create a new location from georeferenced file
  2. Use the Import file tool, or the r.in.gdal module, to import the GeoTIFF files into GRASS GIS.

See also Automated data import below.

Notes

  1. Most Landsat scenes are delivered in north-is-up orientation, hence the import process is straightforward.
  2. If you get ERROR: Input map is rotated - cannot import., the image must be first rotated externally, 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
    

Hint: Minimal disk space copies

Here's a little trick using r.reclass to rename maps (for example, from L71074092_09220040924_B10 to L71074092_09220040924_B1) 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


Automated data import

The following Python script imports Landsat imagery into GRASS' data base. Specifically, the script

  • creates an independent Mapset for each Landsat scene
  • imports and renames bands of a scene as B<id>, e.g. B10, B20, ..., B80.
  • additionaly sets up the timestamp based on MTL metadata file


Note, the (newest) official naming pattern for Landsat scenes -- explained in USGS' Landsat Data Dictionary as the Landsat Scene Identifier -- and all individual bands that compose a scene -- bands have a suffix which is like _B10, _B20, _B30, etc. -- differ from what some Landsat specific modules expect. For example, the modules i.landsat.toar and i.landsat.acca expect the bands to follow a naming pattern such as "scenename.1, .2, .3", etc.


To use the script save it as import_landsat.py file and make sure it is granted the execution permission.

#!/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 any(x in line for x in ('DATE_ACQUIRED', 'ACQUISITION_DATE')):
                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)
        if ('VCID') in ffile:
            name = "".join((os.path.splitext(file)[0].split('_'))[1::2])
        else:
            name = os.path.splitext(file)[0].split('_')[-1]
        if len(name) == 3 and name[-1] == '0':
            band = int(name[1:2])
        elif len(name) == 3 and name[-1] != '0':
            band = int(name[1:3])
        else:
            band = int(name[-1:])
        grass.message('Importing %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 = 'band %d' % band)
        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

After having collected the Landsat scenes of interest in one directory, the script can be used as follows:

  • ./import_landsat.py → the script will walk through a pool directory that contains unique Landsat scene directories (e.g. three directories named after the official naming pattern: LT51800342011158MOR00 LT51810352009079MTI00 LT51820352009326MTI00) and import all bands of each individual scene in their own Mapset
  • ./import_landsat.py LM41890261983200FFF03 → the scrip will import only bands of the specified Landsat scene directory

Post-Processing

Natural color composites

Creating natural (also known as true-) color composites, can be done by

a) using i.colors.enhance to automatically balance the colors for the Red, Green and Blue bands

or

b) equalizing the colors on each of the Red, Green and Blue bands

    r.colors -e map=band1 color=grey
    

Composites can then be produced with the r.composite module.

Resetting the color tables of all bands back to normal greyscale can be done with a for loop:

BASE=L71074092_09220040924

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

See also:

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

BASE=L71074092_09220040924

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

Calculate Top-of-Atmosphere reflectance and band-6 temperature using the i.landsat.toar module. For details, refer to USGS' Landsat Filename and Metadata Changes dedicated page. For Landsat 8, see also here.

 i.landsat.toar input_prefix=$BASE output_prefix=${BASE}_toar metfile=${BASE}_MTL.txt

Note, the resulting temperature map is in Kelvin:

 # convert to degree Celsius
 r.mapcalc "$BASE.temp_celsius = ${BASE}_toar.6 - 273.15"
 r.info -r $BASE.temp_celsius

Haze removal

Simple haze removal can be done with i.landsat.dehaze. This addons applies a bandwise haze correction using tasscap4 (haze) and linear regression.

During the 2000s, prior acquired Landsat data were reprocessed to LPGS (Level 1 Product Generation System). Seems that with this level of processing, haze effect in raw data was removed and now is sufficient to apply i.colors.enhance to obtain a good image with high contrast.

Atmospheric correction

See Atmospheric correction

Cloud identification

Identify clouds in the image with i.landsat.acca:

i.landsat.acca -f input_prefix=226_62_toar. output=226_62.acca_cloudmask

Mask out the clouds:

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

Download sample data

Preprocessed Landsat-7 data for North Carolina

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

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

The above import efforts are not needed since the data are already in a GRASS location.

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

LANDSAT Image classification

See Image classification

Time series analysis

See Time series

Sources