Overview

Test “SegOptim”, an OBIA algorithm for creating segments, and apply afterwards a supervised classification.

Objectives

Familiarize yourself with the algorithm using a test dataset. Afterwards, using a new given dataset, create the input data and perform an OBIA segmentation and classification.

Data

Dataset 1: test data prepared for running the code. Dataset 2: Sentinel images from practice “Change detection”

  • S2A_MSIL2A_20170620T100031_N0205_R122_T33UUP_20170620T100453

  • S2A_MSIL2A_20170829T100031_N0205_R122_T33UUP_20170829T100026

  • Vector of the AOI

Procedure

1. Install and/or load necessary packages and set your working directory

ipak <- function(pkg){
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    if (length(new.pkg)) 
        install.packages(new.pkg, dependencies = TRUE)
    sapply(pkg, require, character.only = TRUE)
}


packages <- c("rgdal","raster","sp","RStoolbox","RColorBrewer","remotes","sf","stars")
ipak(packages)

install_github("joaofgoncalves/SegOptim")
library(SegOptim)
setwd("your/path/")

2. Define paths of input data 

# Path to raster data used for image segmentation
# In SEGM_FEAT directory
inputSegFeat.path <- "test/data/SEGM_FEAT/SegmFeat_WV2_b532.tif"

# Path to training raster data
# [0] Non-invaded areas [1] Acacia dealbata invaded areas
# In TRAIN_AREAS directory
trainData.path <- "test/data/TRAIN_AREAS/TrainAreas.tif"

# Path to raster data used as classification features
# In CLASSIF_FEAT directory
classificationFeatures.path <- c("test/data/CLASSIF_FEAT/ClassifFeat_WV2_NDIs.tif",
                                 "test/data/CLASSIF_FEAT/ClassifFeat_WV2_SpectralBands.tif")

# Path to Orfeo Toolbox binaries 
otbPath <- "D:\\SW\\OTB-7.2.0-Win64\\bin"

## Output file from OTB image segmentation     
outSegmRst.path <- "test/results3/segmRaster.tif"

# The final output file containing the distribution of the target species
outClassRst.path <- "test/results3/WV2_VilarVeiga_AcaciaDealbata_v1.tif"

Check if your input paths exist.

# Check if the files and folders exist! 
if(!file.exists(inputSegFeat.path))
  print("Could not find data for image segmentation!")

if(!file.exists(trainData.path))
  print("Could not find training data!")

if(!(file.exists(classificationFeatures.path[1]) | file.exists(classificationFeatures.path[2])))
  print("Could not find data for classification!")

if(!dir.exists(otbPath))
  print("Could not find Orfeo Toolbox binaries!")

3. Segmentation

This section will use SegOptim interface functions to run OTB’s Large Scale Mean Shift (LSMS) image segmentation algorithm. Segmentation involves the partitioning of an image into a set of jointly exhaustive and mutually disjoint regions (a.k.a. segments, composed by multiple image pixels), that are internally more homogeneous and similar, compared to adjacent ones. Image segments are then related to geographic objects of interest (e.g., forests, agricultural or urban areas) through some form of object-based classification (the following step in this exercise). Source

## Run the segmentation
outSegmRst <- segmentation_OTB_LSMS(
  # Input raster with features/bands to segment
  inputRstPath = inputSegFeat.path, 
  # Algorithm params
  SpectralRange = 3.1, 
  SpatialRange = 4.5, 
  MinSize = 21,
  # Output
  outputSegmRst = outSegmRst.path,
  verbose = TRUE,
  otbBinPath = otbPath,
  lsms_maxiter = 50)

Check that the output paths were created.

# Check the file paths with outputs
print(outSegmRst)

Visualize the segmentation output.

# Load the segmented raster and plot it
segmRst <- raster(outSegmRst$segm)

x11()
plot(segmRst)

4. Load training data and perfom classification

Train data for SegOptim is simply a raster file with digitized training cases. These are later imported into image segments.

# Train data
trainDataRst <- raster(trainData.path)

Classification features are 2 resters that contain:

  • Surface reflectance for the eight WV2 bands;
  • All combinations (without order) of Normalized Difference Indices.
# Classification features
classificationFeatures <- stack(classificationFeatures.path)
# Change the names for each layer
names(classificationFeatures) <- c(paste("NDI_",1:28,sep=""),paste("SpecBand_",1:8,sep=""))
  • calData - is a data frame object containing calibration data for training and evaluating a classifier algorithm.
calData <- prepareCalData(rstSegm = segmRst,
                          trainData = trainDataRst,
                          rstFeatures = classificationFeatures,
                          thresh = 0.5,
                          funs = "mean",
                          minImgSegm = 30,
                          verbose = TRUE)

# Calibrate/evaluate the RF classifier
classifObj <- calibrateClassifier( calData                    = calData,
                                   classificationMethod       = "RF",
                                   balanceTrainData           = FALSE,
                                   balanceMethod              = "ubOver",
                                   evalMethod                 = "10FCV",
                                   evalMetric                 = "Kappa",
                                   minTrainCases              = 30,
                                   minCasesByClassTrain       = 10,
                                   minCasesByClassTest        = 5,
                                   runFullCalibration         = TRUE)
# Get more evaluation measures
evalMatrix <- evalPerformanceClassifier(classifObj)

print(round(evalMatrix,2))

# Finally, predict the class label for the entire image (i.e., outside the training set)
# and also save the classified image
rstPredSegmRF <- predictSegments(classifierObj = classifObj,
                                 calData = calData,
                                 rstSegm = segmRst,
                                 predictFor = "all",
                                 filename = outClassRst.path)
print(rstPredSegmRF)

# Variable importance
varImpPlot(classifObj$ClassObj$FULL)