Commit 800d8ddf authored by carstennh's avatar carstennh
Browse files

add R-package

parent 509bab6b
Package: HaSa
Title: Autonomous Image Sampling and Probability Mapping
Version: 2.0
Authors@R:
person(given = "Carsten",
family = "Neumann",
role = c("aut", "cre"),
email = "carsten.neumann@gfz-potsdam.de",
comment = "developer at GFZ Potsdam")
Description: Calculates samples and related classifiers for mapping gradual probabilities of image units (habitat types). It enables a stepwise discrimination of habitat types in imagery and provides reference sample locations and models for classification of complex scenes.
License: GPL-3
Imports: rgdal, raster, rgeos, spatialEco, randomForest, e1071, velox, leaflet, leafem (<= 0.0.1), htmlwidgets
Remotes:url::https://cran.r-project.org/src/contrib/Archive/velox/velox_0.2.0.tar.gz
Encoding: UTF-8
LazyData: true
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.1.1
Suggests:
knitr,
rmarkdown
VignetteBuilder: knitr
# Generated by roxygen2: do not edit by hand
export(multi_Class_Sampling)
export(plot_Results)
export(write_Out_Samples)
#' Sample Collection for Habitat Types
#'
#'Writes out a set of samples (SpatialPointsDataFrame) into an ESRI shapefile for a selected habitat type. Each point represents a valid sample location that identifies the selected habitat type.
#'
#' @param inPath file path (character) for results of habitat type sampling and probability mapping (same as outPath from function multi_Class_Sampling)
#' @param step step number (numeric)
#' @param className name (character) of habitat type for which samples should be selected
#'
#' @return ESRI shapefile with name: RefHaSa_className_step.shp
#' 1) Point Shape represents pixel that belong to selected habitat type and can be used as reference for further model building
#'
#'
#' @export
###write out selected samples for step = 6 using
write_Out_Samples <- function (inPath, step, className) {
paste(inPath,"step_",step,"_",className,".tif",sep="")
run1<-get(load(paste(inPath,"Run",step,sep="")))
load(paste(inPath,"threshold_step_",step,sep=""))
dummy<-raster::raster(paste(inPath,"step_",step,"_",className,".tif",sep=""))
thres<-threshold[6]
dummy[dummy < thres]<-NA
dummy[dummy >=thres]<-1
collect<-list()
j<-0
###extract only class samples
for ( i in 1:length(run1@ref_samples)) {if (length(dim(run1@ref_samples[[i]])) != 0)
{if (is.na(run1@switch[i]) == F) { j=j+1;collect[[j]]<-run1@ref_samples[[i]][which(run1@ref_samples[[i]]@data==1),]}else
{j=j+1;collect[[j]]<-run1@ref_samples[[i]][which(run1@ref_samples[[i]]@data==2),]}}}
result<-do.call(rbind,collect)
res<-raster::extract(dummy,result)
if (length(which(is.na(res))) >0) { res<-result[-which(is.na(res)),] }
rgdal::writeOGR(res, layer="result", dsn=paste(inPath,"RefHaSa_",className,"_",step,".shp",sep=""), driver="ESRI Shapefile")
}
# HaSa - HabitatSampler
#
# Copyright (C) 2020 Carsten Neumann (GFZ Potsdam, carsten.neumann@gfz-potsdam.de)
#
# This software was developed within the context of the project
# NaTec - KRH (www.heather-conservation-technology.com) funded
# by the German Federal Ministry of Education and Research BMBF
# (grant number: 01 LC 1602A).
# The BMBF supports this project as research for sustainable development (FONA); www.fona.de.
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
###################################################################################
clip<-function(raster,shape) {
raster::rasterOptions(progress="text")
a1_crop<-raster::crop(raster,shape)
step1<-raster::rasterize(shape,a1_crop)
step1<-raster::reclassify(step1,c(1,200,1))
a1_crop*step1
}
#' Sentinel-2A/B satellite time series stack
#'
#'co-registrated and atmospherically corrected (reflectance) 9-waveband scenes for an entire phenological cylce in 2018
#'
#'
#'The area recorded is a open heathland (dominant species:Calluna vulgaris) under different management status on a former military training ground in NE Germany, trees are mainly pine, birch and oak
#'
#' @source Data: ESA European Space Agency
#' @source Processing: GTS2 German Research Centre for Geosciences GFZ
#'
#' @format A RasterBrick object (package:raster) with 54 layers:
#' \describe{
#' \item{Wavebands:}{"blue", "green", "red", "redEdge1", "redEdge2", "redEdge3", "nearInfrared", "shortwaveInfrared1", "shortwaveInfrared2", }
#' \item{Dates:}{"2018-03-18", "2018-04-17", "2018-05-07", "2018-06-06", "2018-07-16", "2018-07-31", "2018-09-09", "2018-09-19", "2018-10-14", }
#' \item{Pixel Size:}{resampled to 10m, }
#' \item{Region:}{Kyritz-Ruppiner Heide, NE Germany, Brandenburg.}
#' }
"Sentinel_Stack_2018"
#' Habitat type point locations
#'
#'Spatial points of habitat types marked on Sentinel_Stack_2018
#'
#'
#'Example data set for habitat types that can be marked as spatial points in a known area of the input image
#'
#' @source German Research Centre for Geosciences GFZ
#'
#' @format A SpatialPointsDataFrame (package:sp) with 7 points:
#' \describe{
#' \item{Data$class:}{deciduous", "coniferous", "heath_young", "heath_old", "heath_shrub", "bare_ground", "xeric_grass", }
#' \item{Projection:}{"+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0". }
#' }
"Example_Reference_Points"
#' Habitat type spectra
#'
#'Habitat types (rows) and corresponding spectral predictors (columns)
#'
#'
#'Example data set for habitat types that can be extracted from imagery using spatial points or that is generated from a reference databases (spectral libary) with spectral predictors that correspond to the layers of the input image
#'
#' @source German Research Centre for Geosciences GFZ
#'
#' @format A Data.Frame with 7 observations and 54 variables:
#' \describe{
#' \item{rows:}{class.names - deciduous", "coniferous", "heath_young", "heath_old", "heath_shrub", "bare_ground", "xeric_grass", }
#' \item{columns:}{spectral predictors either extracted from spectral library with columns = image layers or extracted from spatial point locations. }
#' }
"Example_Reference_Table"
# HaSa - HabitatSampler
#
# Copyright (C) 2020 Carsten Neumann (GFZ Potsdam, carsten.neumann@gfz-potsdam.de)
#
# This software was developed within the context of the project
# NaTec - KRH (www.heather-conservation-technology.com) funded
# by the German Federal Ministry of Education and Research BMBF
# (grant number: 01 LC 1602A).
# The BMBF supports this project as research for sustainable development (FONA); www.fona.de.
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along
# with this program. If not, see <http://www.gnu.org/licenses/>.
###################################################################################
sample_nb<-function(raster,nb_samples,sample_type,nb_mean,nb_it,buffer,reference,model,area,mtry,last,seed,init.seed) {
###
n_channel<-length(raster::names(raster))
###velox
rID=raster[[1]]
rID[]=1:(nrow(rID)*ncol(rID))
r=raster::stack(rID,raster)
ras.vx <- velox::velox(r)
###
l<-1 ###6. opt=260
model1<-1
pbtn1<-matrix(1,nrow=1,ncol=1)
pbtn2<-matrix(2,nrow=1,ncol=1)
m<-vector("numeric",length=length(nb_samples))
layer<-list()
for (r in nb_samples ) {
############################################################################################################################################################
if (last==T) { reference<-rbind(reference,rep(3000,ncol(reference))) }
n<-nb_it
sample_size<-r
max_samples_per_class<-sample_size*5
if (init.seed == "sample") {seed2<-sample(c(1:1000000),size=nb_mean,replace=F)} else {seed2<-init.seed}
oobe<-matrix(NA,nrow=n,ncol=nb_mean)
models<-list()
points<-list()
dif<-matrix(NA,nrow=nb_mean,ncol=nrow(reference))
channel<-matrix(NA,nrow=nb_mean,ncol=nrow(reference))
switch<-matrix(NA,nrow=nb_mean,ncol=nrow(reference))
pb <- utils::txtProgressBar(min=1, max=nb_mean, style=3)
for ( k in 1:nb_mean) {
for ( j in 1:n ) {
###Vorbereitung Klassifizierung
if ( j == 1 ) {classes<-as.factor(c(1,1))
if (sample_type == "random") {set.seed(seed2[k]); pbt<-raster::sampleRandom(raster, size=sample_size,sp=T)}
if (sample_type == "regular") {pbt<-raster::sampleRegular(raster, size=sample_size,sp=T)}
pbt<-spatialEco::point.in.poly(pbt,area)[,1:n_channel]
f<-which(is.na(pbt@data[1]))
if ( length(f) != 0 ) {pbt<-pbt[-f,]}
set.seed(seed2[k]); classes<-as.factor(sample(c(1:2),size=nrow(pbt),replace=T))
if (length(levels(classes)) <2) {break}
data<-as.data.frame(cbind(classes,pbt@data))
}
########################################################################################################################################
if (model == "rf") {
model1<- randomForest::randomForest(as.factor(classes) ~ ., data = data,mtry=mtry)
if (is.na(mean(model1$err.rate[,1])) == TRUE) {break}
oobe[j,k]<-mean(model1$err.rate[,1])}
###
if (model == "svm") {
model1<- e1071::svm(as.factor(classes) ~ ., data = data)
co<-length(which(as.numeric(as.character(model1$fitted))-as.numeric(as.character(classes)) == 0))
if (co == 0) {break}
oobe[j,k]<-1-(co/length(classes))}
#if ( j > 1) {if (oobe[j,k] < 0.02 || abs(oobe[(j-1),k]-oobe[j,k]) <= 0.011 )
if ( j > 1) {if (oobe[j,k] < 0.02)
{models[[k]]<-model1; points[[k]]<-rbind(pbtn1,pbtn2); break}
if (oobe[(j-1),k] <= oobe[j,k]) {models[[k]]<-model_pre; points[[k]]<-rbind(pbtn1_pre,pbtn2_pre); break}
if (j == n & oobe[j,k] >= 0.02) {models[[k]]<-"NULL"; points[[k]]<-"NULL"; break}
}
model_pre<-model1
pbtn1_pre<-pbtn1
pbtn2_pre<-pbtn2
oobe<<-oobe
########################################################################################################################################
if (model == "rf") {correct<-which(as.numeric(as.character(classes)) - as.numeric(as.character(model1$predicted)) == 0)}
if (model == "svm") {correct<-which(as.numeric(as.character(model1$fitted))-as.numeric(as.character(classes)) == 0)}
########################################################################################################################################
if ( length(which(classes[correct] == 1)) == 0 ) {if ( j ==1) {break}else{pbtn1<-pbtn1}}else
{d1<-correct[which(classes[correct] == 1)]
###generate new samples from only correctly classified samples [label 1]
p1<-pbt@coords[d1,]
pbtn1<-as.data.frame(cbind(classes[d1],matrix(p1,ncol=2)))
sp::coordinates(pbtn1)<-c("V2","V3")
sp::proj4string(pbtn1) <- sp::proj4string(pbt)
poly <- rgeos::gBuffer(spgeom=pbtn1, width=buffer,byid=TRUE)
test<-ras.vx$extract(sp=poly)
for ( i in 1:length(test) ){s1<-dim(test[[i]])[1]; if (s1 <= 5) {test[[i]]<-test[[i]] }else { set.seed(seed); test[[i]]<-test[[i]][sample(c(1:s1),5,replace=F),] }}
for ( i in 1:length(test) ) { if ( i == 1 ) { co<-raster::xyFromCell(raster,test[[i]][,1]) }else {co<-rbind(co,raster::xyFromCell(raster,test[[i]][,1]))} }
pbtn1<-as.data.frame(cbind(rep(1,nrow(co)),co))
sp::coordinates(pbtn1)<-c("x","y")
test1<-as.matrix(do.call(rbind,test)[,-1]); if (ncol(test1)==1) {test1<-t(test1)}; colnames(test1)<-raster::names(raster)
if ( length(which(is.na(test1))) > 0 ) {pbtn1<-pbtn1[complete.cases(test1),];test1<-test1[complete.cases(test1),]}
}
if (class(test1)[1] == "numeric") {test1<-t(matrix(test1))}
if (nrow(test1) == 0) {break}
##############################
##############################
if ( length(which(classes[correct] == 2)) == 0 ) {if ( j ==1) {break}else{pbtn2<-pbtn2}}else
{d2<-correct[which(classes[correct] == 2)]
###generate new samples from only correctly classified samples [label 2]
p2<-pbt@coords[d2,]
pbtn2<-as.data.frame(cbind(classes[d2],matrix(p2,ncol=2)))
sp::coordinates(pbtn2)<-c("V2","V3")
sp::proj4string(pbtn2) <- sp::proj4string(pbt)
poly <- rgeos::gBuffer(spgeom=pbtn2, width=buffer,byid=TRUE)
test<-ras.vx$extract(sp=poly)
for ( i in 1:length(test) ){s1<-dim(test[[i]])[1]; if (s1 <= 5) {test[[i]]<-test[[i]] }else { set.seed(seed); test[[i]]<-test[[i]][sample(c(1:s1),5,replace=F),] }}
for ( i in 1:length(test) ) { if ( i == 1 ) { co<-xyFromCell(raster,test[[i]][,1]) }else {co<-rbind(co,xyFromCell(raster,test[[i]][,1]))} }
pbtn2<-as.data.frame(cbind(rep(2,nrow(co)),co))
sp::coordinates(pbtn2)<-c("x","y")
test2<-as.matrix(do.call(rbind,test)[,-1]); if (ncol(test2)==1) {test2<-t(test2)}; colnames(test2)<-raster::names(raster)
if ( length(which(is.na(test2))) > 0 ) {pbtn2<-pbtn2[complete.cases(test2),];test2<-test2[complete.cases(test2),]}
}
if (class(test2)[1] == "numeric") {test2<-t(matrix(test2))}
if (nrow(test2) == 0) {break}
######################################
###Gleichverteilung samples in Klassen
di<-c(nrow(pbtn1),nrow(pbtn2))
if (abs(nrow(pbtn1)-nrow(pbtn2)) > min(di)*0.3) {
if (which.min(di) == 2) {set.seed(seed); d3<-sample(1:nrow(pbtn1),nrow(pbtn2),replace=F); pbtn1<-pbtn1[d3,]; test1<-test1[d3,] }else
{set.seed(seed); d4<-sample(1:nrow(pbtn2),nrow(pbtn1),replace=F);pbtn2<-pbtn2[d4,]; test2<-test2[d4,]}
}
#####################################
###max Klassenbelegungswert
if ( nrow(pbtn1) > max_samples_per_class ) {set.seed(seed); dr<-sample(1:nrow(pbtn1),max_samples_per_class,replace=F)
pbtn1<-pbtn1[dr,]
test1<-test1[dr,]}
if ( nrow(pbtn2) > max_samples_per_class ) {set.seed(seed); dr<-sample(1:nrow(pbtn2),max_samples_per_class,replace=F)
pbtn2<-pbtn2[dr,]
test2<-test2[dr,]
}
########################################################################################################################################
data<-as.data.frame(cbind(append(pbtn1@data$V1,pbtn2@data$V1),rbind(test1,test2))) ##data
names(data)[1]<-"classes"
classes<-data$classes
pbt<-rbind(pbtn1,pbtn2)
}
utils::setTxtProgressBar(pb, k)}
if (length(models) == 0 | length(which(models=="NULL")) == length(models)) {stop("No Models - would you be so kind to increase init.samples, please")}
if ( length(which(models=="NULL")) > 0) {models<-models[-which(models=="NULL")]}
for (jj in 1:nrow(reference)) {ref<-jj;rr=3
for ( i in 1:length(models) ) {
if ( i == 1 ) {dummy<-as.numeric(as.character(stats::predict(models[[i]],newdata=reference)))
if (dummy[ref] != 2) {dummy[dummy==1]<-3;dummy[dummy==2]<-1;dummy[dummy==3]<-2;switch[i,jj]<-i}}else
{dummy2<-as.numeric(as.character(stats::predict(models[[i]],newdata=reference)))
if (dummy2[ref] != 2) {dummy2[dummy2==1]<-3;dummy2[dummy2==2]<-1;dummy2[dummy2==3]<-2;switch[i,jj]<-i}
dummy_set<-dummy
dummy<-dummy+dummy2
dummy3<-dummy/dummy[ref]
dif[i,jj]<-dummy3[ref]-max(dummy3[-ref])
if ( i > 2) {if (dummy3[ref]-max(dummy3[-ref]) > dif[(rr-1),jj]) {dummy<-dummy; channel[i,jj]<-i;dif[(rr-1),jj]<-dif[i,jj]}else {dummy<-dummy_set}}
}}
m[l]<-max(dif[2,],na.rm=T)
}
index<<-which.max(dif[2,])
ch<-as.numeric(na.omit(channel[,index]));if(length(ch) == 0) {stop("No optimal classifier - would you be so kind to adjust init.samples & nb_models, please")}
acc<<-(round(m[l]^2,2)/0.25)
print(paste("class=",index," difference=",(round(m[l]^2,2)/0.25),sep=""))
l<-l+1
}
close(pb)
mod_all<-models
models<-models[ch];print(paste("n_models =",length(models)))
switch<-switch[ch,index]
points<-points[ch]
dif<-dif[2,]
##########################################################################################################################################
###Vohersage
ch<-c(1:length(models))
switch<-switch
j<-1
for ( i in ch ) {
if ( j == 1 ) {result1<-raster::predict(object=raster, model=models[[i]])
if (is.na(switch[i]) == F) {result1<-raster::reclassify(result1,rbind(c(0.5,1.5,2),c(1.6,2.5,1)))}}else
{result1<-raster::stack(result1,raster::predict(object=raster, model=models[[i]]))
if (is.na(switch[i]) == F) {result1[[j]]<-raster::reclassify(result1[[j]],rbind(c(0.5,1.5,2),c(1.6,2.5,1)))}
}
print(j)
j<-j+1}
###
dummy<-raster::brick(result1)
dummy<-raster::calc(dummy,fun=sum)
layer[[1]]<-dummy
setClass("Habitat",representation(models="list", ref_samples="list", switch="vector", layer="list", mod_all="list", class_ind="numeric", seeds="numeric"))
new("Habitat", models = models, ref_samples = points, switch = switch, layer = layer, mod_all = mod_all, class_ind = dif, seeds = seed2)
}
#' Perform Habitat Sampling and Probability Mapping
#'
#'This is the main function that performs everything: specify the input imagery, select model type, initiate sampling and model building, generates interactive maps and produce final probability raster output
#'
#' @param in.raster satellite time series stack (rasterBrickObject) or just any type of image (*rasterObject)
#' @param init.samples starting number of spatial locations
#' @param sample_type distribution of spatial locations c("random","regular")
#' @param nb_models number of models (independent classifiers) to collect
#' @param nb_it number of iterations for model accuracy
#' @param buffer distance (in m) for new sample collection around initial samples (depends on pixel size)
#' @param reference reference spectra either SpatialPointsDataFrame (shape file) or data.frame with lines = classes, column = predictors]
#' @param model which machine learning classifier to use c("rf", "svm") for random forest or suppurt vector machine implementation
#' @param mtry number of predictor used at random forest splitting nodes (mtry << n predictors)
#' @param last only true for one class classifier c("FALSE", TRUE")
#' @param seed set seed for reproducable results
#' @param init.seed "sample" for new or use run1@seeds to reproduce previous steps
#' @param outPath output path for saving results
#' @param step at which step should the procedure start, e.g. use step = 2 if the first habitat is already extracted
#' @param classNames vector with class names in the order of reference spectra
#' @param n_classes total number of classes (habitat types) to be separated
#' @param multiTest number of test runs to compare different probability outputs
#' @param RGB rgb channel numbers for image plot
#'
#' @return 4 files per step:
#' 1) Habitat type probability map as geocoded *.kml layer and *.tif raster files and *.png image output
#' 2) A Habitat object consisting of 7 slots: \cr
#' run1@models - list of selcted classifiers \cr
#' run1@ref_samples - list of SpatialPointsDataFrames with same length as run1@models holding reference labels [1,2] for each selected model \cr
#' run1@switch - vector of lenght run1@models indicating if target class equals 2, if not NA the labels need to be switched \cr
#' run1@layer - raster map of habitat type probability \cr
#' run1@mod_all - list of all classifiers (equals nb_models) \cr
#' run1@class_ind - vector of predictive distance measure for all habitats \cr
#' run1@seeds - vector of seeds for random sampling \cr
#' all files are saved with step number, the *.tif file is additionally saved with class names
#'
#' @examples next steps start automatically, after command line input of:
#' 1) number of the apropriate map if multiTest > 1
#' 2) probability threshold for habitat type extraction
#' 3) decision to sample again y/n
#' 4) adjust starting number of samples and number of models
#'
#' for threshold evaluation an interactive map is plotted in the web browser
#'
#' if convergence fails / no models can be selected / init.samples are to little / or another error occurs, restart next step with:
#' in.raster = out.raster
#' reference = out.reference
#' step = specify next step number
#' classNames = out.names
#'
#' @export
###################################################################################
multi_Class_Sampling<-function(in.raster, init.samples=30, sample_type="regular", nb_models=200, nb_it=10, buffer, reference, model="rf", mtry=10, last=F, seed=3, init.seed="sample", outPath, step=1, classNames, n_classes, multiTest=1,RGB=c(1,2,3)) {
###first steps: data preparation
if (class(reference) == "SpatialPointsDataFrame") { reference<-as.data.frame(raster::extract(in.raster,reference)) }
area <- as(extent(in.raster), 'SpatialPolygons')
area <- sp::SpatialPolygonsDataFrame(area, data.frame( ID=1:length(p)))
sp::proj4string(area)<-sp::proj4string(in.raster)
r<-RGB[1]; g<-RGB[2]; b<-RGB[3]
col<-colorRampPalette(c("lightgrey","orange","yellow","limegreen","forestgreen"))
source("inner_procedure.r")
source("clip,r")
source("plot_interactive.r")
#############################################################################################
r<-n_classes
if (names(in.raster)[1] != colnames(reference)[1]) { colnames(reference)<-names(in.raster) }
if (step != 1) {if (step<11) {load(paste(outPath,paste("threshold_step_0",step-1,sep=""),sep=""))}else{load(paste(outPath,paste("threshold_step_",step-1,sep=""),sep=""))}}
for ( i in step:r) {print(paste(paste("init.samples = ",init.samples),paste("models = ",nb_models))); if ( i == r) {last=T}
if ( multiTest > 1 ) { test<-list(); maFo<-list(); new.names<-list(); decision = "0"
####################################################################################################################
while (decision == "0") {
for ( rs in 1:multiTest ) {
########################
maFo_rf<-sample_nb(raster=in.raster,nb_samples=seq(init.samples,init.samples,init.samples),sample_type=sample_type,nb_mean=nb_models,nb_it=nb_it,buffer=buffer,reference=reference,model=model,area=area,mtry=mtry,last=last,seed=seed,init.seed=init.seed)
########################
maFo[[rs]]<-maFo_rf
test[[rs]]<-maFo_rf@layer[[1]]
new.names[[rs]]<-index
if ( rs == multiTest ) {par(mar=c(2,2,2,3),mfrow=n2mfrow(multiTest))
for ( rr in 1:length(test) ) { plot(test[[rr]], col = col(200),main="", legend.shrink=1); mtext(side=3, paste(rr,classNames[new.names[[rr]]], sep=" "),font =2) }}
}
decision<-readline("Which distribution is acceptable/ or sample again [../0]: ")
}
maFo_rf<-maFo[[as.numeric(decision)]]
index<-new.names[[as.numeric(decision)]]
####################################################################################################################
}else{
########################
maFo_rf<-sample_nb(raster=in.raster,nb_samples=seq(init.samples,init.samples,init.samples),sample_type=sample_type,nb_mean=nb_models,nb_it=nb_it,buffer=buffer,reference=reference,model=model,area=area,mtry=mtry,last=last,seed=seed,init.seed=init.seed)
########################
}
dummy<-maFo_rf@layer[[1]]
iplot(x=dummy,y=a1,HaTy=classNames[index])
decision<-readline("Threshold for Habitat Extraction or Sample Again [../0]: ")
sample2<-init.samples
models2<-nb_models
while (decision == "0") { decision2<-readline("Adjust init.samples/nb.models or auto [../.. or 0]: ")
if (decision2 != "0") {sample2<-as.numeric(strsplit(decision2, split="/")[[1]][1])
models2<-as.numeric(strsplit(decision2, split="/")[[1]][2])}else
{sample2<-sample2+50
models2<-models2+15}
print(paste(paste("init.samples = ",sample2),paste("models = ",models2)))
maFo_rf<-sample_nb(raster=in.raster,nb_samples=seq(sample2,sample2,sample2),sample_type=sample_type,nb_mean=models2,nb_it=nb_it,buffer=buffer,reference=reference,model=model,area=area,mtry=mtry,last=last,seed=seed,init.seed=init.seed)
########################
dummy<-maFo_rf@layer[[1]]
iplot(x=dummy,y=a1,HaTy=classNames[index])
decision<-readline("Threshold for Habitat Extraction or Sample Again [../0]: ")
}
run1<-maFo_rf
if (i <10) {ni<-paste("0",i,sep="")}else{ni<-i}
save(run1,file=paste(outPath,paste("Run",ni,sep=""),sep=""))
raster :: writeRaster(dummy,filename=paste(outPath,paste("step_",ni,paste("_",classNames[index],sep=""),".tif",sep=""),sep=""), format="GTiff")
#savePlot("step_1.png",type="png")
kml<-raster :: projectRaster(dummy, crs="+proj=longlat +datum=WGS84", method='ngb')
raster :: KML(kml,paste(outPath,paste("step_",ni,sep=""),sep=""))
thres<-as.numeric(decision)
dummy<-maFo_rf@layer[[1]]
dummy[dummy < thres]<-1
dummy[dummy >=thres]<-NA
reference<-reference[-index,]; out.reference<<-reference
classNames<-classNames[-index]; out.names<<-classNames
in.raster<-in.raster*dummy; out.raster<<-in.raster
print(paste(paste("Habitat",i),"Done"))
if( i == r ) {print("Congratulation - you finally made it towards the last habitat"); break()}
colnames(reference)<-names(in.raster)
if ( i == 1) {threshold<-thres; save(threshold,file=paste(outPath,paste("threshold_step_",ni,sep=""),sep=""))}else{threshold<-append(threshold,thres); save(threshold,file=paste(outPath,paste("threshold_step_",ni,sep=""),sep=""))}
decision2<-readline("Adjust init.samples/nb.models or auto [../.. or 0]: ")
if (decision2 != "0") {init.samples<-as.numeric(strsplit(decision2, split="/")[[1]][1])
nb_models<-as.numeric(strsplit(decision2, split="/")[[1]][2])}else {init.samples<-init.samples+50
nb_models<-nb_models+15}
}
}
# HaSa - HabitatSampler
#
# Copyright (C) 2020 Carsten Neumann (GFZ Potsdam, carsten.neumann@gfz-potsdam.de)
#
# This software was developed within the context of the project
# NaTec - KRH (www.heather-conservation-technology.com) funded
# by the German Federal Ministry of Education and Research BMBF
# (grant number: 01 LC 1602A).
# The BMBF supports this project as research for sustainable development (FONA); www.fona.de.
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.