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.
We will use two Landsat 8 (OLI) scenes clipped to the nc sample data region
Download the shared data and let us create a new mapset in the "nc_spm_08_grass7" location to do further analysis.
First we will create a new mapset called "user1_L8" .To make use of the existing data in the NC location, we will add the existing 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
# Launch GRASS GIS, -c creates new mapset user1_l8
# If in Windows, use Power shell
grass72 D:/grassdata/nc_spm_08_grass7/user1_L8/ -c
# If in Windows Power shell, type in bash and enter
bash
# 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.
# Change directory to the input Landsat 8 data
cd /d/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"
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.
# 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
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.
# 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"
Now create false colour and true colour composite for better visualization
# 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"
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.
# 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
Task: Visually compare the cloud coverage with the false color composite
We will compute NDVI and NDWI from the spectral bands using the map calculator
.
# 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"
Now classify the LC80150352016200LGN00 (cloud free) using unsupervised sequential Maxlike algorithm.
Main steps are:
# 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.
Last changed: 02 Nov 2017
GRASS GIS manual main index | Topics index | Keywords Index | Full index | Raster index | Vector index | Temporal index |
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.