LULCC Analysis in NNP
Yujue Wang
2021/1/10
LULCC Analysis of Namadgi National Park
Operations carried out in this research in R include: 1) Landsat images pre-processing 2) Train the random forest classifier 3) Test model 4) Classify LULC using fitted model
Part 1: Landsat images preprocessing
Adopted NDVI, NDWI, MNDWI and MSAVI indices instead of spectral bands for image classification.
NDVI - Normalised Difference Vegetation Index = (nir - red)/(nir + red)
NDWI - Normalised Difference Water Index = (green - nir)/(green + nir)
MNDWI - Modified Normalised Difference Water Index = (green - swir2)/(green + swir2)
MSAVI - Modified Soil Adjusted Vegetation Index = nir + 0.5 - (0.5 * sqrt((2 * nir + 1)^2 - 8 * (nir - (2 * red))))
Landsat images and DEM data were downloaded in USGS website: https://glovis.usgs.gov/
The shapefile of study area was downloaded in ACTmapi website:http://www.actmapi.act.gov.au/download.html
#clear environment and console
rm(list = ls())
cat("\014")
#Load required packages
library(raster)
library(rgdal)
library(caret)
library(RStoolbox)# import data
beginCluster()## 8 cores detected, using 7meta2003 <- readMeta("rawdata/LT05_L1TP_090085_20031117_20161204_01_T1/LT05_L1TP_090085_20031117_20161204_01_T1_MTL.txt")
meta2013 <- readMeta("rawdata/LC08_L1TP_090085_20131230_20170427_01_T1/LC08_L1TP_090085_20131230_20170427_01_T1_MTL.txt")
meta2019 <- readMeta("rawdata/LC08_L1TP_090085_20190113_20190131_01_T1/LC08_L1TP_090085_20190113_20190131_01_T1_MTL.txt")
meta2020 <- readMeta("rawdata/LC08_L1TP_090085_20201201_20201217_01_T1/LC08_L1TP_090085_20201201_20201217_01_T1_MTL.txt")
metatrain <- readMeta("rawdata/LC08_L1TP_090085_20201014_20201104_01_T1/LC08_L1TP_090085_20201014_20201104_01_T1_MTL.txt")
summary(meta2003)## Scene: LT50900852003321ASA03
## Satellite: LANDSAT5
## Sensor: TM
## Date: 2003-11-17
## Path/Row: 90/85
## Projection: +proj=utm +zone=55 +datum=WGS84 +units=m +no_defs
##
## Data:
## FILES QUANTITY CATEGORY
## B1_dn LT05_L1TP_090085_20031117_20161204_01_T1_B1.TIF dn image
## B2_dn LT05_L1TP_090085_20031117_20161204_01_T1_B2.TIF dn image
## B3_dn LT05_L1TP_090085_20031117_20161204_01_T1_B3.TIF dn image
## B4_dn LT05_L1TP_090085_20031117_20161204_01_T1_B4.TIF dn image
## B5_dn LT05_L1TP_090085_20031117_20161204_01_T1_B5.TIF dn image
## B6_dn LT05_L1TP_090085_20031117_20161204_01_T1_B6.TIF dn image
## B7_dn LT05_L1TP_090085_20031117_20161204_01_T1_B7.TIF dn image
## QA_dn LT05_L1TP_090085_20031117_20161204_01_T1_BQA.TIF dn qa
##
## Available calibration parameters (gain and offset):
## dn -> radiance (toa)
## dn -> reflectance (toa)
## dn -> brightness temperature (toa)summary(meta2013)## Scene: LC80900852013364LGN01
## Satellite: LANDSAT8
## Sensor: OLI_TIRS
## Date: 2013-12-30
## Path/Row: 90/85
## Projection: +proj=utm +zone=55 +datum=WGS84 +units=m +no_defs
##
## Data:
## FILES QUANTITY CATEGORY
## B1_dn LC08_L1TP_090085_20131230_20170427_01_T1_B1.TIF dn image
## B2_dn LC08_L1TP_090085_20131230_20170427_01_T1_B2.TIF dn image
## B3_dn LC08_L1TP_090085_20131230_20170427_01_T1_B3.TIF dn image
## B4_dn LC08_L1TP_090085_20131230_20170427_01_T1_B4.TIF dn image
## B5_dn LC08_L1TP_090085_20131230_20170427_01_T1_B5.TIF dn image
## B6_dn LC08_L1TP_090085_20131230_20170427_01_T1_B6.TIF dn image
## B7_dn LC08_L1TP_090085_20131230_20170427_01_T1_B7.TIF dn image
## B9_dn LC08_L1TP_090085_20131230_20170427_01_T1_B9.TIF dn image
## B10_dn LC08_L1TP_090085_20131230_20170427_01_T1_B10.TIF dn image
## B11_dn LC08_L1TP_090085_20131230_20170427_01_T1_B11.TIF dn image
## B8_dn LC08_L1TP_090085_20131230_20170427_01_T1_B8.TIF dn pan
## QA_dn LC08_L1TP_090085_20131230_20170427_01_T1_BQA.TIF dn qa
##
## Available calibration parameters (gain and offset):
## dn -> radiance (toa)
## dn -> reflectance (toa)
## dn -> brightness temperature (toa)summary(meta2019)## Scene: LC80900852019013LGN00
## Satellite: LANDSAT8
## Sensor: OLI_TIRS
## Date: 2019-01-13
## Path/Row: 90/85
## Projection: +proj=utm +zone=55 +datum=WGS84 +units=m +no_defs
##
## Data:
## FILES QUANTITY CATEGORY
## B1_dn LC08_L1TP_090085_20190113_20190131_01_T1_B1.TIF dn image
## B2_dn LC08_L1TP_090085_20190113_20190131_01_T1_B2.TIF dn image
## B3_dn LC08_L1TP_090085_20190113_20190131_01_T1_B3.TIF dn image
## B4_dn LC08_L1TP_090085_20190113_20190131_01_T1_B4.TIF dn image
## B5_dn LC08_L1TP_090085_20190113_20190131_01_T1_B5.TIF dn image
## B6_dn LC08_L1TP_090085_20190113_20190131_01_T1_B6.TIF dn image
## B7_dn LC08_L1TP_090085_20190113_20190131_01_T1_B7.TIF dn image
## B9_dn LC08_L1TP_090085_20190113_20190131_01_T1_B9.TIF dn image
## B10_dn LC08_L1TP_090085_20190113_20190131_01_T1_B10.TIF dn image
## B11_dn LC08_L1TP_090085_20190113_20190131_01_T1_B11.TIF dn image
## B8_dn LC08_L1TP_090085_20190113_20190131_01_T1_B8.TIF dn pan
## QA_dn LC08_L1TP_090085_20190113_20190131_01_T1_BQA.TIF dn qa
##
## Available calibration parameters (gain and offset):
## dn -> radiance (toa)
## dn -> reflectance (toa)
## dn -> brightness temperature (toa)summary(meta2020)## Scene: LC80900852020336LGN00
## Satellite: LANDSAT8
## Sensor: OLI_TIRS
## Date: 2020-12-01
## Path/Row: 90/85
## Projection: +proj=utm +zone=55 +datum=WGS84 +units=m +no_defs
##
## Data:
## FILES QUANTITY CATEGORY
## B1_dn LC08_L1TP_090085_20201201_20201217_01_T1_B1.TIF dn image
## B2_dn LC08_L1TP_090085_20201201_20201217_01_T1_B2.TIF dn image
## B3_dn LC08_L1TP_090085_20201201_20201217_01_T1_B3.TIF dn image
## B4_dn LC08_L1TP_090085_20201201_20201217_01_T1_B4.TIF dn image
## B5_dn LC08_L1TP_090085_20201201_20201217_01_T1_B5.TIF dn image
## B6_dn LC08_L1TP_090085_20201201_20201217_01_T1_B6.TIF dn image
## B7_dn LC08_L1TP_090085_20201201_20201217_01_T1_B7.TIF dn image
## B9_dn LC08_L1TP_090085_20201201_20201217_01_T1_B9.TIF dn image
## B10_dn LC08_L1TP_090085_20201201_20201217_01_T1_B10.TIF dn image
## B11_dn LC08_L1TP_090085_20201201_20201217_01_T1_B11.TIF dn image
## B8_dn LC08_L1TP_090085_20201201_20201217_01_T1_B8.TIF dn pan
## QA_dn LC08_L1TP_090085_20201201_20201217_01_T1_BQA.TIF dn qa
##
## Available calibration parameters (gain and offset):
## dn -> radiance (toa)
## dn -> reflectance (toa)
## dn -> brightness temperature (toa)summary(metatrain)## Scene: LC80900852020288LGN00
## Satellite: LANDSAT8
## Sensor: OLI_TIRS
## Date: 2020-10-14
## Path/Row: 90/85
## Projection: +proj=utm +zone=55 +datum=WGS84 +units=m +no_defs
##
## Data:
## FILES QUANTITY CATEGORY
## B1_dn LC08_L1TP_090085_20201014_20201104_01_T1_B1.TIF dn image
## B2_dn LC08_L1TP_090085_20201014_20201104_01_T1_B2.TIF dn image
## B3_dn LC08_L1TP_090085_20201014_20201104_01_T1_B3.TIF dn image
## B4_dn LC08_L1TP_090085_20201014_20201104_01_T1_B4.TIF dn image
## B5_dn LC08_L1TP_090085_20201014_20201104_01_T1_B5.TIF dn image
## B6_dn LC08_L1TP_090085_20201014_20201104_01_T1_B6.TIF dn image
## B7_dn LC08_L1TP_090085_20201014_20201104_01_T1_B7.TIF dn image
## B9_dn LC08_L1TP_090085_20201014_20201104_01_T1_B9.TIF dn image
## B10_dn LC08_L1TP_090085_20201014_20201104_01_T1_B10.TIF dn image
## B11_dn LC08_L1TP_090085_20201014_20201104_01_T1_B11.TIF dn image
## B8_dn LC08_L1TP_090085_20201014_20201104_01_T1_B8.TIF dn pan
## QA_dn LC08_L1TP_090085_20201014_20201104_01_T1_BQA.TIF dn qa
##
## Available calibration parameters (gain and offset):
## dn -> radiance (toa)
## dn -> reflectance (toa)
## dn -> brightness temperature (toa)# stack bands along with the metadata
TM2003 <- stackMeta(meta2003)
OLI2013 <- stackMeta(meta2013)
OLI2019 <- stackMeta(meta2019)
OLI2020 <- stackMeta(meta2020)
OLItrain <- stackMeta(metatrain)
# import mask of study area
mask <- readOGR("rawdata/Namadgi National Park Border.shp")
# convert projection
mask <- spTransform(mask, CRS(proj4string(TM2003)))
# import DEM
dem <- raster ("rawdata/DEM.tif")
# convert projection
dem <- projectRaster(dem,TM2003)
# clip to the area mask
dem_clip1 <- crop(dem, extent(mask))
dem_clip <- mask(dem_clip1, mask)
TM2003_clip1 <- crop(TM2003, extent(mask))
TM2003_clip <- mask(TM2003_clip1, mask)
OLI2013_clip1 <- crop(OLI2013, extent(mask))
OLI2013_clip <- mask(OLI2013_clip1, mask)
OLI2019_clip1 <- crop(OLI2019, extent(mask))
OLI2019_clip <- mask(OLI2019_clip1, mask)
OLI2020_clip1 <- crop(OLI2020, extent(mask))
OLI2020_clip <- mask(OLI2020_clip1, mask)
OLItrain_clip1 <- crop(OLItrain, extent(mask))
OLItrain_clip <- mask(OLItrain_clip1, mask)
# Atmospheric correction using DOS
#(correct DN to at-surface-reflectance with DOS (Chavez decay model))
TM2003_dos <- radCor(TM2003_clip, meta2003, method = "dos", darkProp = 0.01)
OLI2013_dos <- radCor(OLI2013_clip, meta2013, method = "dos", darkProp = 0.01)
OLI2019_dos <- radCor(OLI2019_clip, meta2019, method = "dos", darkProp = 0.01)
OLI2020_dos <- radCor(OLI2020_clip, meta2020, method = "dos", darkProp = 0.01)
OLItrain_dos <- radCor(OLItrain_clip, metatrain, method = "dos", darkProp = 0.01)
# Topographic correction
TM2003_dos_topcor <-topCor(TM2003_dos, dem_clip, metaData = meta2003, method = "C")
OLI2013_dos_topcor <- topCor(OLI2013_dos, dem_clip, metaData = meta2013, method = "C")
OLI2019_dos_topcor <- topCor(OLI2019_dos, dem_clip, metaData = meta2019, method = "C")
OLI2020_dos_topcor <- topCor(OLI2020_dos, dem_clip, metaData = meta2020, method = "C")
OLItrain_dos_topcor <- topCor(OLItrain_dos, dem_clip, metaData = metatrain, method = "C")SpectralIndices function can calculate spectral indices automatically by entering the corresponding band number.
Note that the band numbers of Landsat 5 TM and Landsat 8 OLI are different.
# Calculate Spectral Indices
SI2003 <- spectralIndices(TM2003_dos_topcor, blue = 1, green = 2, red = 3, nir = 4, swir2 = 5, swir3 = 7)
SI2013 <- spectralIndices(OLI2013_dos_topcor, blue = 2, green = 3, red = 4, nir = 5,swir2 = 6,swir3 = 7)
SI2019 <- spectralIndices(OLI2019_dos_topcor, blue = 2, green = 3, red = 4, nir = 5,swir2 = 6,swir3 = 7)
SI2020 <- spectralIndices(OLI2020_dos_topcor, blue = 2, green = 3, red = 4, nir = 5,swir2 = 6,swir3 = 7)
SItrain <- spectralIndices(OLItrain_dos_topcor, blue = 2, green = 3, red = 4, nir = 5,swir2 = 6,swir3 = 7)
# SI2020@data@names
# Output Spectral Indices
writeRaster(SI2003[[7]], "classifyLULC/MNDWI_2003.TIF", overwrite=TRUE)
writeRaster(SI2003[[8]], "classifyLULC/MSAVI_2003.TIF", overwrite=TRUE)
writeRaster(SI2003[[11]], "classifyLULC/NDVI_2003.TIF", overwrite=TRUE)
writeRaster(SI2003[[13]], "classifyLULC/NDWI_2003.TIF", overwrite=TRUE)
writeRaster(SI2013[[7]], "classifyLULC/MNDWI_2013.TIF", overwrite=TRUE)
writeRaster(SI2013[[8]], "classifyLULC/MSAVI_2013.TIF", overwrite=TRUE)
writeRaster(SI2013[[11]], "classifyLULC/NDVI_2013.TIF", overwrite=TRUE)
writeRaster(SI2013[[13]], "classifyLULC/NDWI_2013.TIF", overwrite=TRUE)
writeRaster(SI2019[[7]], "classifyLULC/MNDWI_2019.TIF", overwrite=TRUE)
writeRaster(SI2019[[8]], "classifyLULC/MSAVI_2019.TIF", overwrite=TRUE)
writeRaster(SI2019[[11]], "classifyLULC/NDVI_2019.TIF", overwrite=TRUE)
writeRaster(SI2019[[13]], "classifyLULC/NDWI_2019.TIF", overwrite=TRUE)
writeRaster(SI2020[[7]], "classifyLULC/MNDWI_2020.TIF", overwrite=TRUE)
writeRaster(SI2020[[8]], "classifyLULC/MSAVI_2020.TIF", overwrite=TRUE)
writeRaster(SI2020[[11]], "classifyLULC/NDVI_2020.TIF", overwrite=TRUE)
writeRaster(SI2020[[13]], "classifyLULC/NDWI_2020.TIF", overwrite=TRUE)
writeRaster(SItrain[[7]], "classifyLULC/MNDWI_train.TIF", overwrite=TRUE)
writeRaster(SItrain[[8]], "classifyLULC/MSAVI_train.TIF", overwrite=TRUE)
writeRaster(SItrain[[11]], "classifyLULC/NDVI_train.TIF", overwrite=TRUE)
writeRaster(SItrain[[13]], "classifyLULC/NDWI_train.TIF", overwrite=TRUE)
### Output pre-processed raster
writeRaster(TM2003_dos_topcor, "classifyLULC/preprocessed_TM2003_dos_topcor.TIF", format="GTiff")
writeRaster(OLI2013_dos_topcor, "classifyLULC/preprocessed_OLI2013_dos_topcor.TIF", format="GTiff")
writeRaster(OLI2019_dos_topcor, "classifyLULC/preprocessed_OLI2019_dos_topcor.TIF", format="GTiff")
writeRaster(OLI2020_dos_topcor, "classifyLULC/preprocessed_OLI2020_dos_topcor.TIF", format="GTiff")
writeRaster(OLItrain_dos_topcor, "classifyLULC/preprocessed_OLItrain_dos_topcor.TIF", format="GTiff")
endCluster()Part 2: Train the random forest classifier
The main task is to fit random forest model with training data
#import training data
ndvi <- raster("classifyLULC/NDVI_train.tif")
ndwi <- raster("classifyLULC/NDWI_train.tif")
mndwi <- raster("classifyLULC/MNDWI_train.tif")
msavi <- raster("classifyLULC/MSAVI_train.tif")
img <- stack(ndvi, ndwi, mndwi, msavi)
names(img) <- c ("NDVI", "NDWI", "MNDWI", "MSAVI")The training image is converted from original Landsat image after pre-processing and supervised classification into the shapefile format.
(using Semi-Automatic Classification Plugin version 6 in QGIS and referring to the 2010 land use image from GLOBELAND30 website: http://www.globallandcover.com/home.html)
#import classified shapefile
trainData_0 <- shapefile("rawdata/traindata.shp")
trainData <- spTransform(trainData_0, CRS(proj4string(ndvi)))
responseCol <- "class"
# match class with img
dfAll = data.frame(matrix(vector(), nrow = 0, ncol = length(names(img)) + 1))
for (i in 1:length(unique(trainData[[responseCol]]))){
category <- unique(trainData[[responseCol]])[i]
categorymap <- trainData[trainData[[responseCol]] == category,]
dataSet <- extract(img, categorymap)
if(is(trainData, "SpatialPointsDataFrame")){
dataSet <- cbind(dataSet, class = as.numeric(rep(category, nrow(dataSet))))
dfAll <- rbind(dfAll, dataSet[complete.cases(dataSet),])
}
if(is(trainData, "SpatialPolygonsDataFrame")){
dataSet <- dataSet[!unlist(lapply(dataSet, is.null))]
dataSet <- lapply(dataSet, function(x){cbind(x, class = as.numeric(rep(category, nrow(x))))})
df <- do.call("rbind", dataSet)
dfAll <- rbind(dfAll, df)
}
}
#write.table(dfAll, file = "matchclass_train", )
#Fit model
nsamples <- 20000
sdfAll <- dfAll[sample(1:nrow(dfAll), nsamples), ]
modFit_rf <- train(as.factor(class) ~ NDVI+NDWI+MNDWI+MSAVI, method = "rf", data = sdfAll, na.action = na.exclude)
beginCluster()
preds_rf_train <- clusterR(img, raster::predict, args = list(model = modFit_rf))
endCluster()
output = preds_rf_train
plot(preds_rf_train)
writeRaster(preds_rf_train, "classifyLULC/predictedRFclassificationtrain.TIF", overwrite=TRUE )
Part 3: Test model
Overall accuracy and kappa coefficient are calculated in QGIS using Grass plugin.
# calculate kappa coefficient (doesn't work)
#library(sp)
#train_raster <- raster("rawdata/traindata.tif")
#s = 11
#r.kappa <- raster.change(preds_rf_train, train_raster, d = s, stat = "kappa", mask = TRUE)
#plot(r.kappa, main="Kappa")Part 4: Classify LULC using fitted model
classify LULC in 2003, 2013, 2019 and 2020 based on fitted model
# import spectral indices of 2003
ndvi2003 <- raster("classifyLULC/NDVI_2003.tif")
ndwi2003 <- raster("classifyLULC/NDWI_2003.tif")
mndwi2003 <- raster("classifyLULC/MNDWI_2003.tif")
msavi2003 <- raster("classifyLULC/MSAVI_2003.tif")
img2003 <- stack(ndvi2003, ndwi2003, mndwi2003, msavi2003)
names(img2003) <- c ("NDVI", "NDWI", "MNDWI", "MSAVI")
# classify LULC of 2003
beginCluster()
preds_rf_2003 <- clusterR(img2003, raster::predict, args = list(model = modFit_rf))
endCluster()
output = preds_rf_2003
plot(preds_rf_2003)
writeRaster(preds_rf_2003, "classifyLULC/predictedRFclassification2003.TIF", overwrite=TRUE )
The rest operations are similar
# import spectral indices of 2013
ndvi2013 <- raster("classifyLULC/NDVI_2013.tif")
ndwi2013 <- raster("classifyLULC/NDWI_2013.tif")
mndwi2013 <- raster("classifyLULC/MNDWI_2013.tif")
msavi2013 <- raster("classifyLULC/MSAVI_2013.tif")
img2013 <- stack(ndvi2013, ndwi2013, mndwi2013, msavi2013)
names(img2013) <- c ("NDVI", "NDWI", "MNDWI", "MSAVI")
# classify LULC of 2013
beginCluster()
preds_rf_2013 <- clusterR(img2013, raster::predict, args = list(model = modFit_rf))
endCluster()
output = preds_rf_2013
plot(preds_rf_2013)
writeRaster(preds_rf_2013, "classifyLULC/predictedRFclassification2013.TIF", overwrite=TRUE )
# import spectral indices of 2019
ndvi2019 <- raster("classifyLULC/NDVI_2019.tif")
ndwi2019 <- raster("classifyLULC/NDWI_2019.tif")
mndwi2019 <- raster("classifyLULC/MNDWI_2019.tif")
msavi2019 <- raster("classifyLULC/MSAVI_2019.tif")
img2019 <- stack(ndvi2019, ndwi2019, mndwi2019, msavi2019)
names(img2019) <- c ("NDVI", "NDWI", "MNDWI", "MSAVI")
# classify LULC of 2019
beginCluster()
preds_rf_2019 <- clusterR(img2019, raster::predict, args = list(model = modFit_rf))
endCluster()
output = preds_rf_2019
plot(preds_rf_2019)
writeRaster(preds_rf_2019, "classifyLULC/predictedRFclassification2019.TIF", overwrite=TRUE )
# import spectral indices of 2020
ndvi2020 <- raster("classifyLULC/NDVI_2020.tif")
ndwi2020 <- raster("classifyLULC/NDWI_2020.tif")
mndwi2020 <- raster("classifyLULC/MNDWI_2020.tif")
msavi2020 <- raster("classifyLULC/MSAVI_2020.tif")
img2020 <- stack(ndvi2020, ndwi2020, mndwi2020, msavi2020)
names(img2020) <- c ("NDVI", "NDWI", "MNDWI", "MSAVI")
# classify LULC of 2020
beginCluster()
preds_rf_2020 <- clusterR(img2020, raster::predict, args = list(model = modFit_rf))
endCluster()
output = preds_rf_2020
plot(preds_rf_2020)
writeRaster(preds_rf_2020, "classifyLULC/predictedRFclassification2020.TIF", overwrite=TRUE )
Using CA-Markov model to simulate future scenarios
Step 1: Convert raster images to ESRI ASCii format and import them into IDRISI Selva software.
example
Step 2: Calculate Markov probability transition matrix based on LULC images of 2003 and 2013.
example
Step 3: Calculate area transition matrix based on LULC image of 2003 and 2020.
Step 4: Making suitability atlas of 2003 and 2020.
Since there are less constraint factors in natural systems, no additional restrictions are made.
(the scope of water body remains unchanged).
example
Step 5: Prediction LULC in 2019 (2003 as the starting point) using CA-Markov model.
example
Step 6: Test model validation.
Calculate overall accuracy and kappa coefficient of simulated 2019 LULC and the real image (same as before).
Step 7: Predict LULC in 2025, 2035 and 2050 using CA-Markov model and export.
example