Remote sensing analysis in GRASS GIS

display composites

Basic overview of this session

In this session, we will go through the pre-processing of Landsat 8 images covering the topics from importing to deriving vegetation indices and unsupervised classification.

Course content

Data

We will use two Landsat 8 (OLI) scenes clipped to the nc sample data region

Spectral range
Spectral bands of Landsat 8 OLI, (Source: https://landsat.gsfc.nasa.gov/landsat-data-continuity-mission/ )

Data download

For this exercise, download the clipped landsat 8 scenes from here

Workflow overview

How do we work with satellite data in GRASS GIS?

Assuming we do not have our maps in the GRASSDBASE (i.e., grassdata/location/mapset), the first step is to import the datasets to GRASS GIS.

Here is a list of major GRASS GIS modules we will be using in this exercise and the links to their manuals.

  1. Start page of the GRASS GIS manual
  2. Raster processing modules
  3. Vector processing modules
  4. Imagery processing modules
  5. GRASS GIS addons
  6. General commands to start GRASS GIS, set region, list data: grass72, g.region, g.list
  7. Import raster, vector data, import and project on the fly: r.in.gdal, v.in.ogr, r.import, v.import
  8. Raster map info, map calculator, univariate statistics, aggregation: r.info, r.mapcalc, r.univar, r.series,
  9. Vector map info, data base select, univariate statistics, vector to raster: v.info, v.db.select, v.univar, v.to.rast,

Remote sensing analysis in GRASS GIS

Download the shared data and let us create a new mapset in the "nc_spm_08_grass7" location to do further analysis.

To make use of the existing data in the NC location, we will add the mapset "landsat" in the search path. We can run the commands from the command line or use the main GUI and copy the corresponding commands for future replication of the workflow. Note the "Copy" button in the GUI of each module. Maybe you might want to get help about the options and flags of the different commands, e.g. r.colors --help

        
grass72 $HOME/grassdata/nc_spm_08_grass7/user1_l8/ -c
g.proj -p
g.mapsets -p
g.mapsets mapset=landsat operation=add
g.mapsets -p
g.list rast
g.region rast=lsat7_2002_20 res=30 -a
        
        
# Launch GRASS GIS, -c creates new mapset user1_l8
grass72 $HOME/grassdata/nc_spm_08_grass7/user1_l8/ -c
# Let us check the projection of the location
g.proj -p
# List all the mapsets in the search path
g.mapsets -p
# Add the mapset landsat to the search path
g.mapsets mapset=landsat operation=add
# List all the mapsets in the search path
g.mapsets -p
# List all the raster maps in all the mapsets in the search path
g.list type=rast
# Set the computational region 
g.region rast=lsat7_2002_20 res=30 -a
        
    

We will now import the landsat 8 raw data to the newly created mapset.

        
cd $HOME/data_dir/LC80150352016168LGN00
BASE="LC80150352016168LGN00"
for i in "1" "2" "3" "4" "5" "6" "7" "9" "QA" "10" "11"; do
	r.import input=${BASE}_B${i}.TIF output=${BASE}_B${i} resolution=value resolution_value=30
done
r.import input=${BASE}_B8.TIF output=${BASE}_B8 resolution=value resolution_value=15
        
        
# Change directory to the input Landsat 8 data
cd $HOME/data_dir/LC80150352016168LGN00
# Define a variable
BASE="LC80150352016168LGN00"
# Define a loop to import all the bands
for i in "1" "2" "3" "4" "5" "6" "7" "9" "QA" "10" "11"; do
	r.import input=${BASE}_B${i}.TIF output=${BASE}_B${i} resolution=value resolution_value=30
done
# PAN band 8 imported separately because of different spatial resolution
r.import input=${BASE}_B8.TIF output=${BASE}_B8 resolution=value resolution_value=15
        
    

Task: Note that we are using r.import instead of r.in.gdal to import the data. Check the difference between two and explain why we used r.import here?

Task: Repeat the import step for the second scene "LC80150352016200LGN00"

Digital Number (DN) to Top Of Atmosphere (TOA) reflectance

The next step is to convert the digital number (Landsat 8 OLI sensor provides 16 bit data with range between 0 and 65536) to TOA reflectance. For the thermal bands 10 and 11, DN is converted to TOA Brightness Temperature. In GRASS GIS i.landsat.toar can do this step for all the landsat sensors.

        
i.landsat.toar input=${BASE}_B output=${BASE}_toar metfile=${BASE}_MTL.txt sensor=oli8
g.list rast map=. pattern=${BASE}_toar*
        
        
# Convert from DN to TOA reflectance and Brightness Temperature
i.landsat.toar input=${BASE}_B output=${BASE}_toar metfile=${BASE}_MTL.txt sensor=oli8
g.list rast map=. pattern=${BASE}_toar*
        
   

Task: Repeat the conversion step for the second scene "LC80150352016200LGN00"

Task: Set the color palette of Band 10 of LC80150352016200LGN00 to "kelvin" and visulaize with map features

display ST
TOA Temperature from band 10 of Landsat 8 dated 18 July 2016

Data fusion/Pansharpening

Now let us use the PAN band 8 (15 m resolution) to downsample other spectral bands to 15 m resolution. We use an addon i.fusion.hpf which applies a high pass filter addition method to down sample. Here we introduce the long list of addons in GRASS GIS and demonstrate how to install and use them. Check g.extension to install the addons and GRASS GIS addons for the list of available addons.

        
g.region rast=lsat7_2002_20 res=15 -a
g.extension extension=i.fusion.hpf op=add
i.fusion.hpf -l -c pan=${BASE}_toar8 msx=${BASE}_toar1,${BASE}_toar2,${BASE}_toar3,${BASE}_toar4,${BASE}_toar5,${BASE}_toar6,${BASE}_toar7 center=high modulation=max trim=0.0 --o
g.list rast map=. pattern=${BASE}_toar*.hpf
        
        
# Set the region
g.region rast=lsat7_2002_20 res=15 -a
# Install the reqquired addon
g.extension extension=i.fusion.hpf op=add
# Apply the fusion based on high pass filter
i.fusion.hpf -l -c pan=${BASE}_toar8 msx=${BASE}_toar1,${BASE}_toar2,${BASE}_toar3,${BASE}_toar4,${BASE}_toar5,${BASE}_toar6,${BASE}_toar7 center=high modulation=max trim=0.0 --o
# list the fused maps
g.list rast map=. pattern=${BASE}_toar*.hpf
        

Task: Repeat the fusion step for the second scene "LC80150352016200LGN00"

Image Composites

Now create false colour and true colour composite for better visualization

        
g.region rast=lsat7_2002_20 res=15 -a
i.colors.enhance red="${BASE}_toar4.hpf" green="${BASE}_toar3.hpf" blue="${BASE}_toar2.hpf" strength=95
r.composite red="${BASE}_toar4.hpf" green="${BASE}_toar3.hpf" blue="${BASE}_toar2.hpf" output="${BASE}_toar.hpf_comp_432"
i.colors.enhance red="${BASE}_toar5.hpf" green="${BASE}_toar4.hpf" blue="${BASE}_toar3.hpf" strength=95
r.composite red="${BASE}_toar5.hpf" green="${BASE}_toar4.hpf" blue="${BASE}_toar3.hpf" output="${BASE}_toar.hpf_comp_543"
        
        
# Set the region
g.region rast=lsat7_2002_20 res=15 -a
# Enhance the colors in the clipped region
i.colors.enhance red="${BASE}_toar4.hpf" green="${BASE}_toar3.hpf" blue="${BASE}_toar2.hpf" strength=95
# Create RGB composites
r.composite red="${BASE}_toar4.hpf" green="${BASE}_toar3.hpf" blue="${BASE}_toar2.hpf" output="${BASE}_toar.hpf_comp_432"
# Enhance the colors in the clipped region
i.colors.enhance red="${BASE}_toar5.hpf" green="${BASE}_toar4.hpf" blue="${BASE}_toar3.hpf" strength=95
# Create RGB composites
r.composite red="${BASE}_toar5.hpf" green="${BASE}_toar4.hpf" blue="${BASE}_toar3.hpf" output="${BASE}_toar.hpf_comp_543"
        
 

Task: Create the composites for the second scene "LC80150352016200LGN00"

display composite
True color and False color composites of the Landsat 8 image dated 18 July 2016

Cloud mask from the QA layer

Landsat 8 provides quality layer which contains 16bit integer values that represent "bit-packed combinations of surface, atmosphere, and sensor conditions that can affect the overall usefulness of a given pixel". We can use the addon i.landsat8.qc to develop masks. More information about Landsat 8 quality band is given here.

        
g.region rast=lsat7_2002_20 res=15 -a
g.extension extension=i.landsat8.qc op=add
i.landsat8.qc cloud="Maybe,Yes" output=Cloud_Mask_rules.txt
r.reclass input=${BASE}_BQA output=${BASE}_Cloud_Mask rules=Cloud_Mask_rules.txt
r.report -e map=${BASE}_Cloud_Mask units=k -n
        
        
# Set the region
g.region rast=lsat7_2002_20 res=15 -a
# Install the required extension
g.extension extension=i.landsat8.qc op=add
# Create a rule set
i.landsat8.qc cloud="Maybe,Yes" output=Cloud_Mask_rules.txt
# Reclass the BQA band based on the rule set created 
r.reclass input=${BASE}_BQA output=${BASE}_Cloud_Mask rules=Cloud_Mask_rules.txt
# Report the area covered by Cloud
r.report -e map=${BASE}_Cloud_Mask units=k -n
        
    

display cloud mask
False color composite and the derived cloud mask of the Landsat 8 image dated 16 June 2016

Task: Visually compare the cloud coverage with the false color composite

Vegetation Indices

We will compute NDVI and NDWI from the spectral bands using the map calculator

.
        
g.region rast=lsat7_2002_20 res=15 -a
r.mask rast=${BASE}_Cloud_Mask
r.mapcalc "${BASE}_NDVI = (${BASE}_toar5.hpf - ${BASE}_toar4.hpf) / (${BASE}_toar5.hpf + ${BASE}_toar4.hpf) * 1.0"
r.colors ${BASE}_NDVI color=ndvi
r.mapcalc "${BASE}_NDWI = (${BASE}_toar5.hpf - ${BASE}_toar6.hpf) / (${BASE}_toar5.hpf + ${BASE}_toar6.hpf) * 1.0"
r.colors ${BASE}_NDWI color=ndwi
r.mask -r
        
        
# Set the region
g.region rast=lsat7_2002_20 res=15 -a
# Set the cloud mask to avoid computing over clouds
r.mask rast=${BASE}_Cloud_Mask
# Compute NDVI
r.mapcalc "${BASE}_NDVI = (${BASE}_toar5.hpf - ${BASE}_toar4.hpf) / (${BASE}_toar5.hpf + ${BASE}_toar4.hpf) * 1.0"
# Set the color palette
r.colors ${BASE}_NDVI color=ndvi
# Compute NDWI
r.mapcalc "${BASE}_NDWI = (${BASE}_toar5.hpf - ${BASE}_toar6.hpf) / (${BASE}_toar5.hpf + ${BASE}_toar6.hpf) * 1.0"
# Set the color palette
r.colors ${BASE}_NDWI color=ndwi
# Remove the mask
r.mask -r
        
 

Task: Compute the vegetation indices from the second scene "LC80150352016200LGN00"

display composite
False color composite and the derived cloud mask of the Landsat 8 image dated 16 June 2016

Unsupervised Classification

Now classify the LC80150352016200LGN00 (cloud free) using unsupervised sequential Maxlike algorithm.

Main steps are:

        
g.list rast map=. pattern=${BASE}_toar*.hpf
i.group group=${BASE}_hpf subgroup=${BASE}_hpf input=`g.list rast map=. pattern=${BASE}_toar*.hpf sep=","`
i.cluster group=${BASE}_hpf subgroup=${BASE}_hpf sig=${BASE}_hpf classes=8 separation=0.5
i.maxlik group=${BASE}_hpf subgroup=${BASE}_hpf sig=${BASE}_hpf output=${BASE}_hpf.class rej=${BASE}_hpf.rej
        
        
# List the bands needed for classification
g.list rast map=. pattern=${BASE}_toar*.hpf
# add maps to an imagery group for easier management
i.group group=${BASE}_hpf subgroup=${BASE}_hpf input=`g.list rast map=. pattern=${BASE}_toar*.hpf sep=","`
# statistics for unsupervised classification
i.cluster group=${BASE}_hpf subgroup=${BASE}_hpf sig=${BASE}_hpf classes=8 separation=0.5
# Maximum Likelihood unsupervised classification
i.maxlik group=${BASE}_hpf subgroup=${BASE}_hpf sig=${BASE}_hpf output=${BASE}_hpf.class rej=${BASE}_hpf.rej
        
 

Check the manual of i.cluster and i.maxlik.

Other (very) useful links

References


Last changed: 02 Nov 2017

GRASS GIS manual main index | Topics index | Keywords Index | Full index | Raster index | Vector index | Temporal index |

Creative Commons License
Hands-on to GIS and RS with GRASS GIS by Veronica Andreo, Sajid Pareeth and Paulo van Breugel is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License - Thanks to Vaclav Petras for the style.