####################PREMIERES ETAPES################################################
# ---Clear the workspace 
rm(list=ls())
 
#Load libraries
library(raster)
library(rgdal)
library(stringr)
 
###################DEFINIR DES ENTREES###############################################
#set working directory. Il faut definir ici le directoire. Il faut mettre / au lieu de \
wd <- "C:/Users/olofsson/Desktop/STR_example/"
setwd(wd)
 
#Set the project code
Proj <- "MAK2"
# read a raster, GeoTiff or something. Il faut mettre le nom du raster de stratification. 
# Il faut que le raster soit dans le working directory
Raster <- raster('Col_strat_buffer_subset.tif')
 
##Resolution de la carte utilise
Res <- 30
 
#File name with number of samples per stratum. Fichier avec une colonne "map_value" 

rp <- read.csv('STR.csv', header = T)
 
# Fichier de sortie avec les echantillons
output_file <- "STR_selection.csv"
 
# The ISO country code. Code ISO du pays
countrycode <- 'MDG'
 
##Nom du SHP avec les points qu'on va creer
Shp <- "Points_echantillonnage"
 
###NE TOUCHER CE QUI EST ECRIT CI DESSUS######################################################
###NE TOUCHER CE QUI EST ECRIT CI DESSUS######################################################
############################CALCUL DU NOMBRE DE PIXELS PER CLASSE DE LA CARTE##################
t <- Raster@data@attributes
map_area <- as.data.frame(t)
map_area$Ha
map_area$Ha <- map_area$Count*Res*Res/10000
 
#We paste the map_area in the rp
rp$map_area <- map_area$Count
 
 
##########################SELECTION DES ECHANTILLONS#########################################3
#beginCluster()
 
############### Generate 10x times the number of points from overall sample
 
rand_sample <- data.frame(sampleRandom(Raster,(sum(rp$final) *10 + log((sum(rp$map_area)))),xy=TRUE))
names(rand_sample) <- c("x_coord","y_coord","map_value")
              
 
rand_sample$id     <- row(rand_sample)[,1]
rp2 <- merge(rp,data.frame(table(rand_sample$map_value)),by.x="map_value",by.y="Var1",all.x=T)  
rp2[is.na(rp2)]<-0
    
############### Create the list of classes that need to be specifically sampled
to_rtp <- rp2[rp2$Freq <  rp2$final,]$map_value
    
############### Create the list of classes that are enough represented in the random sampling
to_spl <- rp2[rp2$Freq >= rp2$final,]$map_value
    
############### Sample points from the first class
i = 1
    
final <- rand_sample[ rand_sample$id %in% 
                        sample(rand_sample[rand_sample$map_value %in% c(to_spl[i],to_rtp[i]),]$id,
                               rp2[rp2$map_value %in% c(to_spl[i],to_rtp[i]),]$final),]
    
############### Loop into the well represented classes, sample and append
if(length(to_spl) > 1){
for(i in 2:length(to_spl)){
tmp <- rand_sample[ rand_sample$id  %in%  
                      sample(rand_sample[rand_sample$map_value == to_spl[i],]$id, rp2[rp2$map_value == to_spl[i],]$final),]
        final <- rbind(final,tmp)
      }
    }
    
############### Loop into the subrepresented classes, raster_to_point then append
if(length(to_rtp) > 0){
for(i in 1:length(to_rtp)){
        tmp_rtp <- as.data.frame(rasterToPoints(Raster,fun=function(rast){rast==to_rtp[i]}))
        names(tmp_rtp) <- c("x_coord","y_coord","map_value")
        tmp_rtp$id<-row(tmp_rtp)[,1]
        sampling <- min(rp2[rp2$map_value == to_rtp[i],]$final,
                        rp2[rp2$map_value == to_rtp[i],]$map_area)
        tmp<-tmp_rtp[tmp_rtp$id 
                     %in% 
                       sample(tmp_rtp[tmp_rtp$map_value == to_rtp[i],]$id,
                              sampling
                       ),
                     ]
        final <- rbind(final,tmp)                              
      }
    }
 
#endCluster()
    
    points <- final
 
 
##################################################################################################################################
############### Create vector layer with the points
 
        sp_df<-SpatialPointsDataFrame(
          coords=points[,c(1,2)],
          data=data.frame(points[,c(3)]),
          proj4string=CRS(proj4string(Raster))
        )
        
        sp_df <- spTransform(sp_df,CRS("+proj=longlat +datum=WGS84"))
 
 
#################################ON VA CREER LE CSV###############################3
coord <- sp_df@coords
coord.sp <- SpatialPoints(coord)
coord.df <- as.data.frame(coord)
coord.spdf <- SpatialPointsDataFrame(coord.sp, coord.df)
 
#download province boundaries for the country
adm <- getData ('GADM', country= countrycode, level=1)
#match the coordinate systems for the sample points and the boundaries
proj4string(coord.spdf) <-proj4string(adm)
adm1 <- over(coord.spdf, adm)
nsamples <- nrow(coord)
 
id <- matrix(sample(1:nsamples , nsamples , replace=F),nrow = nsamples , ncol =1, dimnames= list(NULL,c("ID")))
id <- as.data.frame(id)
id$ID <- paste(Proj,"-",id$ID)
YCoordinate <- coord[,2]
XCoordinate <- coord[,1]
elevation <- getData("alt", country = countrycode)
slope <- terrain(elevation, opt = "slope")
aspect <- terrain(elevation, opt = "aspect")
elevation <- extract(elevation, cbind(coord[,1], coord[,2]))
slope <- extract(slope, cbind(coord[,1], coord[,2]))
aspect <- extract(aspect, cbind(coord[,1], coord[,2]))
map_code <- sp_df@data[, 1]
 
 
#write CSV file, this can be used directly in Collect Earth
m <- data.frame(id, YCoordinate, XCoordinate, elevation, slope, aspect, map_code) #, cadrillage)
 
write.csv(m, file= output_file, row.names= FALSE) 
 
# write it out to a shapefile for further processing
writeOGR(obj = sp_df, dsn = Shp, layer = Shp, driver="ESRI Shapefile") # this is in geographical projection