R Notebook

This note is intended to use for geospatial data in R

Load packages

library(sf); library(raster); library(ggplot2); library(rgdal); library(geojsonsf); library(tidyr); library(dplyr); library(mlr)

Load training data in R

setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Training_Sample")

The working directory was changed to C:/Users/DELL/OneDrive - tuaf.edu.vn/Seasonal_Rice_Cropland_Mapping_Paper/Data/Training_Sample inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

train_st<-st_read(dsn = "C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Landsat\\Training", layer="Training19")

Error: Cannot open "C:\Users\DELL\OneDrive - tuaf.edu.vn\Seasonal_Rice_Cropland_Mapping_Paper\Data\Landsat\Training"; The source could be corrupt or not supported. See `st_drivers()` for a list of supported formats.

Removing NA values and organize training data and create spatial point shapefile

2013

# 2013
setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Rice_training")

LS_2013<-read.csv("LS_2013.csv",header = T)

LS_2013<-na.omit(LS_2013)

dim(LS_2013); table(LS_2013$Class)

LS_2013$Class[LS_2013$Class=="Cropland"]<-"Rice"

ggplot(LS_2013)+geom_boxplot(aes(x=Class,y=NDVI,fill=Class))

table(LS_2013$Class)

Urban<-sample_n(LS_2013[LS_2013$Class=="Urban",],1300) 

Forest<-sample_n(LS_2013[LS_2013$Class=="Forest",],1300) 

LS_2013$Class[LS_2013$Class=="Urban"]<-NA

LS_2013$Class[LS_2013$Class=="Forest"]<-NA

LS_2013<-na.omit(LS_2013)

LS_2013<-rbind(LS_2013,Urban, Forest); table(LS_2013$Class)

LS_2013$Class[LS_2013$Class=="Forest"]<-"NonRice"
LS_2013$Class[LS_2013$Class=="Urban"]<-"NonRice"
LS_2013$Class[LS_2013$Class=="Water"]<-"NonRice"
table(LS_2013$Class)
# Create point shapefile

LS_2013<-st_as_sf(cbind.data.frame(LS_2013,geojsonsf::geojson_sf(LS_2013$.geo)))

st_write(LS_2013,dsn= "Training_2013.gpkg",driver = "ESRI Shapefile", append = T)

2014

# 2014
setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Rice_training")

LS_2014<-read.csv("LS_2014.csv",header = T)

LS_2014$Class[LS_2014$Class=="Cropland"]<-"Rice"

LS_2014<-na.omit(LS_2014)

ggplot(LS_2014)+geom_boxplot(aes(x=Class,y=NDVI,fill=Class)); table(LS_2014$Class)

Urban<-sample_n(LS_2014[LS_2014$Class=="Urban",],1200) # Sample 2100 points from urban feature 
Forest<-sample_n(LS_2014[LS_2014$Class=="Forest",],1300) 

LS_2014$Class[LS_2014$Class=="Urban"]<-NA

LS_2014$Class[LS_2014$Class=="Forest"]<-NA

LS_2014<-na.omit(LS_2014)

LS_2014<-rbind(LS_2014,Urban, Forest); table(LS_2014$Class)

LS_2014$Class[LS_2014$Class=="Forest"]<-"NonRice"
LS_2014$Class[LS_2014$Class=="Urban"]<-"NonRice"
LS_2014$Class[LS_2014$Class=="Water"]<-"NonRice"
table(LS_2014$Class)

# Create point shapefile

LS_2014<-st_as_sf(cbind.data.frame(LS_2014,geojsonsf::geojson_sf(LS_2014$.geo)))

st_write(LS_2014,dsn= "Training_2014.gpkg",driver = "ESRI Shapefile", append = T)

##2015

# 2015
setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Rice_training")

LS_2015<-read.csv("LS_2015.csv",header = T)

LS_2015<-na.omit(LS_2015)

LS_2015$Class[LS_2015$Class=="Cropland"]<-"Rice"

ggplot(LS_2015)+geom_boxplot(aes(x=Class,y=NDVI,fill=Class))

Urban<-sample_n(LS_2015[LS_2015$Class=="Urban",],1200) # Sample 2100 points from urban feature 

Forest<-sample_n(LS_2015[LS_2015$Class=="Forest",],1300) 

LS_2015$Class[LS_2015$Class=="Urban"]<-NA

LS_2015$Class[LS_2015$Class=="Forest"]<-NA

LS_2015<-na.omit(LS_2015)

LS_2015<-rbind(LS_2015,Urban, Forest); table(LS_2015$Class)

LS_2015$Class[LS_2015$Class=="Forest"]<-"NonRice"
LS_2015$Class[LS_2015$Class=="Urban"]<-"NonRice"
LS_2015$Class[LS_2015$Class=="Water"]<-"NonRice"
table(LS_2015$Class)

# Create point shapefile

LS_2015<-st_as_sf(cbind.data.frame(LS_2015,geojsonsf::geojson_sf(LS_2015$.geo)))

st_write(LS_2015,dsn= "Training_2015.gpkg",driver = "ESRI Shapefile", append = T)

##2016

# 2016
setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Rice_training")
LS_2016<-read.csv("LS_2016.csv",header = T)

LS_2016<-na.omit(LS_2016)

LS_2016$Class[LS_2016$Class=="Cropland"]<-"Rice"

ggplot(LS_2016)+geom_boxplot(aes(x=Class,y=NDVI,fill=Class))

Urban<-sample_n(LS_2016[LS_2016$Class=="Urban",],1300) # Sample 2100 points from urban feature 

Forest<-sample_n(LS_2016[LS_2016$Class=="Forest",],1300) 

LS_2016$Class[LS_2016$Class=="Urban"]<-NA

LS_2016$Class[LS_2016$Class=="Forest"]<-NA

LS_2016<-na.omit(LS_2016)

LS_2016<-rbind(LS_2016,Urban, Forest); table(LS_2016$Class)

LS_2016$Class[LS_2016$Class=="Forest"]<-"NonRice"
LS_2016$Class[LS_2016$Class=="Urban"]<-"NonRice"
LS_2016$Class[LS_2016$Class=="Water"]<-"NonRice"
table(LS_2016$Class)
# Create point shapefile

LS_2016<-st_as_sf(cbind.data.frame(LS_2016,geojsonsf::geojson_sf(LS_2016$.geo)))

st_write(LS_2016,dsn= "Training_2016.gpkg",driver = "ESRI Shapefile", append = T)

2017

# 2017
setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Rice_training")
LS_2017<-read.csv("LS_2017.csv",header = T)

LS_2017<-na.omit(LS_2017)

LS_2017$Class[LS_2017$Class=="Cropland"]<-"Rice"

ggplot(LS_2017)+geom_boxplot(aes(x=Class,y=NDVI,fill=Class))

Urban<-sample_n(LS_2017[LS_2017$Class=="Urban",],1200) # Sample 2100 points from urban feature 

Forest<-sample_n(LS_2017[LS_2017$Class=="Forest",],1200) 

LS_2017$Class[LS_2017$Class=="Urban"]<-NA

LS_2017$Class[LS_2017$Class=="Forest"]<-NA

LS_2017<-na.omit(LS_2017)

LS_2017<-rbind(LS_2017,Urban, Forest); table(LS_2017$Class)

LS_2017$Class[LS_2017$Class=="Forest"]<-"NonRice"
LS_2017$Class[LS_2017$Class=="Urban"]<-"NonRice"
LS_2017$Class[LS_2017$Class=="Water"]<-"NonRice"
table(LS_2017$Class)

# Create point shapefile

LS_2017<-st_as_sf(cbind.data.frame(LS_2017,geojsonsf::geojson_sf(LS_2017$.geo)))

st_write(LS_2017,dsn= "Training_2017.gpkg",driver = "ESRI Shapefile", append = T)

2018

# 2018
setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Rice_training")
LS_2018<-read.csv("LS_2018.csv",header = T)

LS_2018<-na.omit(LS_2018)

LS_2018$Class[LS_2018$Class=="Cropland"]<-"Rice"

ggplot(LS_2018)+geom_boxplot(aes(x=Class,y=NDVI,fill=Class))

Urban<-sample_n(LS_2018[LS_2018$Class=="Urban",],1200) # Sample 2100 points from urban feature 

Forest<-sample_n(LS_2018[LS_2018$Class=="Forest",],1100) 

LS_2018$Class[LS_2018$Class=="Urban"]<-NA

LS_2018$Class[LS_2018$Class=="Forest"]<-NA

LS_2018<-na.omit(LS_2018)

LS_2018<-rbind(LS_2018,Urban, Forest)

LS_2018$Class[LS_2018$Class=="Forest"]<-"NonRice"
LS_2018$Class[LS_2018$Class=="Urban"]<-"NonRice"
LS_2018$Class[LS_2018$Class=="Water"]<-"NonRice"
table(LS_2018$Class)
# Create point shapefile

LS_2018<-st_as_sf(cbind.data.frame(LS_2018,geojsonsf::geojson_sf(LS_2018$.geo)))

st_write(LS_2018,dsn= "Training_2018.gpkg",driver = "ESRI Shapefile", append = T)

2019

# 2019
setwd("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Rice_training")
LS_2019<-read.csv("LS_2019.csv",header = T)

LS_2019<-na.omit(LS_2019)

LS_2019$Class[LS_2019$Class=="Cropland"]<-"Rice"

ggplot(LS_2019)+geom_boxplot(aes(x=Class,y=NDVI,fill=Class))

Urban<-sample_n(LS_2019[LS_2019$Class=="Urban",],1400) # Sample 2100 points from urban feature 

Forest<-sample_n(LS_2019[LS_2019$Class=="Forest",],1093) 

LS_2019$Class[LS_2019$Class=="Urban"]<-NA

LS_2019$Class[LS_2019$Class=="Forest"]<-NA

LS_2019<-na.omit(LS_2019)

LS_2019<-rbind(LS_2019,Urban, Forest); table(LS_2019$Class)

LS_2019$Class[LS_2019$Class=="Forest"]<-"NonRice"
LS_2019$Class[LS_2019$Class=="Urban"]<-"NonRice"
LS_2019$Class[LS_2019$Class=="Water"]<-"NonRice"
table(LS_2019$Class)
# Create point shapefile

LS_2019<-st_as_sf(cbind.data.frame(LS_2019,geojsonsf::geojson_sf(LS_2019$.geo)))

st_write(LS_2019,dsn= "Training_2019.gpkg",driver = "ESRI Shapefile", append = T)
                                    

Random Forest with mlr package

2013

library(mlr); library(randomForest); library(parallel); library(parallelMap)
# Preprocessing data 
RF_2013<-LS_2013 %>% data.frame() %>% select(1:9,-6,16,17,24)

RF_2013$Class<-as.factor(RF_2013$Class)
# Create a task and learner 
LS_Task <- makeClassifTask(data = RF_2013, target = "Class")

LS_learner<-makeLearner("classif.randomForest")

# Search for hyperparameter turnings 

forestParamSpace <- makeParamSet(
  makeIntegerParam("ntree", lower = 100, upper = 500),
  makeIntegerParam("mtry", lower = 4, upper = 12),
  makeIntegerParam("nodesize", lower = 1, upper =10),
  makeIntegerParam("maxnodes", lower = 5, upper = 30)) # Set a range of values to search for.

randSearch <- makeTuneControlRandom(maxit = 100)

cvForTuning <- makeResampleDesc("CV", iters = 5)

parallelStartSocket(cpus = detectCores())

tunedForestPars <- tuneParams(LS_learner, task = LS_Task,
                            resampling = cvForTuning,
                            par.set = forestParamSpace,
                            control = randSearch)
parallelStop()

tunedForestPars
# Traing the random forest
tunedForest <- setHyperPars(LS_learner, par.vals = tunedForestPars$x)

tunedForestModel <- train(tunedForest, LS_Task)

# Cross-validate the model using the hyperparameter turnings 
## outer <- makeResampleDesc("CV", iters = 5)

## forestWrapper <- makeTuneWrapper("classif.randomForest", resampling = cvForTuning, par.set = forestParamSpace, control = randSearch)

## parallelStartSocket(cpus = detectCores())

## cvWithTuning <- resample(forestWrapper, LS_Task, resampling = outer)

## parallelStop()

## cvWithTuning
# Cross-validating the model without hyperparameter turnings 

# kFold <- makeResampleDesc(method = "RepCV", folds = 10, reps = 50, stratify = TRUE)

# kFoldCV <- resample(learner = LS_learner, task = LS_Task, resampling = kFold, measures = list(mmce, acc))
# Check the confusion matrix 

# calculateConfusionMatrix(kFoldCV$pred, relative = TRUE)

Tune result: Op. pars: ntree=182; mtry=9; nodesize=2; maxnodes=30 mmce.test.mean=0.1148977

Try

# Preprocessing data 
library(caTools)

RF_2013<-LS_2018 %>% data.frame() %>% select(1:30,-6)

RF_2013$Class<-as.factor(RF_2013$Class)

# Split data in two sets 
split<-sample.split(RF_2013$Class, SplitRatio = 0.8)

# training set
training_set<-subset(RF_2013, split==T)

test_set<-subset(RF_2013, split==F)

table(test_set$Class)
# Create a task and learner 
RF_Task <- makeClassifTask(data = training_set, target = "Class")

RF_learner<-makeLearner("classif.randomForest")

RF_train<-train(RF_learner,RF_Task) 

# Holdout cross-validation
set.seed(234)
RF_10_folds<-makeResampleDesc(method="RepCV", stratify = T, folds=10, reps=5)

RF_cross_validation <- resample(learner = RF_learner, task = RF_Task, resampling = RF_10_folds, measures = list(mmce, acc)) 
# Cross-validation 

calculateConfusionMatrix(RF_cross_validation$pred,relative = T)

# Independent testing
test_set$Class<-as.factor(test_set$Class)

LS_pred1<- predict(RF_train, newdata=test_set[,-1],type="class")

LS<-data.frame(Class=LS_pred1$data)

names(LS)<-"Class"

library(e1071)
library(caret)
test_set$Class<-as.factor(test_set$Class)

cm2<-confusionMatrix(LS$Class,test_set$Class)
cm2
library(raster)

library(sp)

library(rgeos)

library(rgdal)

Calculate the areas

2013

library(raster)
# provinces shapefiles

folder<-"C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Areas"

provinces<-st_read(dsn=folder,layer = "Red_River_Delta")

pro<-st_transform(provinces,crs = 32648)

# Read raster 
LS<-raster("C:\\Users\\DELL\\OneDrive - tuaf.edu.vn\\Seasonal_Rice_Cropland_Mapping_Paper\\Data\\Landsat\\Classified_Nomasked\\map2019.tif")

NAvalue(LS)<-0

names(LS)<-"Classified_2013"
# Clip raster by pprovinces

myarea<-function(image,vectors){
  for (i in 1:nrow(vectors)){
    if (vectors$VARNAME_1[i]=="Bac Ninh"){
      BN<-vectors[i,2]
      BN_crop<-crop(image,BN,snap="out")
      BN_mask<-mask(BN_crop,BN)
      df1<-as.data.frame(BN_mask) %>%   group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(BN_mask)[1] * res(BN_mask)[2]/10000)
      df1<-df1[1,]
      df1$Province<-"Bac Ninh"
      print(df1)

    } else if(vectors$VARNAME_1[i]=="Ha Nam"){
      Hanam<-vectors[i,2]
      Hanam_crop<-crop(image,Hanam,snap="out")
      Hanam_mask<-mask(Hanam_crop,Hanam)
      df2<-as.data.frame(Hanam_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(Hanam_mask)[1] * res(Hanam_mask)[2]/10000)
      df2<-df2[1,]
      df2$Province<-"Ha Nam"
      print(df2)
    }else if(vectors$VARNAME_1[i]=="Ha Noi"){
      HN<-vectors[i,2]
      HN_crop<-crop(image,HN,snap="out")
      HN_mask<-mask(HN_crop,HN)
      df3<-as.data.frame(HN_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(HN_mask)[1] * res(HN_mask)[2]/10000)
      df3<-df3[1,]
      df3$Province<-"Ha Noi"
      print(df3)
    }else if(vectors$VARNAME_1[i]=="Hai Duong"){
      HD<-vectors[i,2]
      HD_crop<-crop(image,HD,snap="out")
      HD_mask<-mask(HD_crop,HD)
      df4<-as.data.frame(HD_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(HD_mask)[1] * res(HD_mask)[2]/10000)
      df4<-df4[1,]
      df4$Province<-"Hai Duong"
      print(df4)
    } else if (vectors$VARNAME_1[i]=="Hai Phong"){
      HP<-vectors[i,2]
      HP_crop<-crop(image,HP,snap="out")
      HP_mask<-mask(HP_crop,HP)
      df5<-as.data.frame(HP_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(HP_mask)[1] * res(HP_mask)[2]/10000)
      df5<-df5[1,]
      df5$Province<-"Hai Phong"
      print(df5)
    }else if(vectors$VARNAME_1[i]=="Hung Yen"){
      HY<-vectors[i,2]
      HY_crop<-crop(image,HY,snap="out")
      HY_mask<-mask(HY_crop,HY)
      df6<-as.data.frame(HY_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(HY_mask)[1] * res(HY_mask)[2]/10000)
      df6<-df6[1,]
      df6$Province<-"Hung Yen"
      print(df6)
    }else if(vectors$VARNAME_1[i]=="Nam Dinh"){
     ND<-vectors[i,2]
     ND_crop<-crop(image,ND,snap="out")
     ND_mask<-mask(ND_crop,ND)
     df7<-as.data.frame(ND_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(ND_mask)[1] * res(ND_mask)[2]/10000) 
     df7<-df7[1,]
     df7$Province<-"Nam Dinh"
     print(df7)
    }else if(vectors$VARNAME_1[i]=="Ninh Binh"){
     NB<-vectors[i,2]
     NB_crop<-crop(image,NB,snap="out")
     NB_mask<-mask(NB_crop,NB)
     df8<-as.data.frame(NB_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(NB_mask)[1] * res(NB_mask)[2]/10000) 
     df8<-df8[1,]
     df8$Province<-"Ninh Binh"
     print(df8)
    }else if(vectors$VARNAME_1[i]=="Thai Binh"){
     TB<-vectors[i,2]
     TB_crop<-crop(image,TB,snap="out")
     TB_mask<-mask(TB_crop,TB)
     df9<-as.data.frame(TB_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(TB_mask)[1] * res(TB_mask)[2]/10000) 
     df9<-df9[1,]
     df9$Province<-"Thai Binh"
     print(df9)
    }  else{
      VP<-vectors[i,2]
      VP_crop<-crop(image,VP,snap="out")
      VP_mask<-mask(VP_crop,VP)
      df10<-as.data.frame(VP_mask) %>% group_by(Classified_2013) %>% tally() %>% mutate(area = n * res(VP_mask)[1] * res(VP_mask)[2]/10000) 
      df10<-df10[1,]
      df10$Province<-"Vinh Phuc"
      print(df10)
    }}} 

df<-myarea(LS,pro)