<jupyter_start><jupyter_text>Define working environment    The following cells are used to: - Import needed libraries- Set the environment variables for Python, Anaconda, GRASS GIS and R statistical computing - Define the ["GRASSDATA" folder](https://grass.osgeo.org/grass73/manuals/helptext.html), the name of "location" and "mapset" where you will to work. **Import libraries**<jupyter_code>## Import libraries needed for setting parameters of operating system 
import os
import sys

## Import library for temporary files creation 
import tempfile 

## Import Pandas library
import pandas as pd

## Import Numpy library
import numpy

## Import Psycopg2 library (interection with postgres database)
import psycopg2 as pg

# Import Math library (usefull for rounding number, e.g.)
import math

## Import Subprocess + subprocess.call
import subprocess
from subprocess import call, Popen, PIPE, STDOUT<jupyter_output><empty_output><jupyter_text>Environment variables when working on Linux Mint    **Set 'Python' and 'GRASS GIS' environment variables** Here, we set [the environment variables allowing to use of GRASS GIS](https://grass.osgeo.org/grass64/manuals/variables.html) inside this Jupyter notebook. Please change the directory path according to your own system configuration.<jupyter_code>### Define GRASS GIS environment variables for LINUX UBUNTU Mint 18.1 (Serena)
# Check is environmental variables exists and create them (empty) if not exists.
if not 'PYTHONPATH' in os.environ:
    os.environ['PYTHONPATH']=''
if not 'LD_LIBRARY_PATH' in os.environ:
    os.environ['LD_LIBRARY_PATH']=''
# Set environmental variables
os.environ['GISBASE'] = '/home/tais/SRC/GRASS/grass_trunk/dist.x86_64-pc-linux-gnu'
os.environ['PATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'bin')
os.environ['PATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'script')
os.environ['PATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'lib')
#os.environ['PATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'etc','python')
os.environ['PYTHONPATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'etc','python')
os.environ['PYTHONPATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'etc','python','grass')
os.environ['PYTHONPATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'etc','python','grass','script')
os.environ['PYTHONLIB'] = '/usr/lib/python2.7'
os.environ['LD_LIBRARY_PATH'] += os.pathsep + os.path.join(os.environ['GISBASE'],'lib')
os.environ['GIS_LOCK'] = '$$'
os.environ['GISRC'] = os.path.join(os.environ['HOME'],'.grass7','rc')
os.environ['PATH'] += os.pathsep + os.path.join(os.environ['HOME'],'.grass7','addons')
os.environ['PATH'] += os.pathsep + os.path.join(os.environ['HOME'],'.grass7','addons','bin')
os.environ['PATH'] += os.pathsep + os.path.join(os.environ['HOME'],'.grass7','addons')
os.environ['PATH'] += os.pathsep + os.path.join(os.environ['HOME'],'.grass7','addons','scripts')

## Define GRASS-Python environment
sys.path.append(os.path.join(os.environ['GISBASE'],'etc','python'))<jupyter_output><empty_output><jupyter_text>**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-** **Display current environment variables of your computer**<jupyter_code>## Display the current defined environment variables
for key in os.environ.keys():
    print "%s = %s \t" % (key,os.environ[key])<jupyter_output>MDMSESSION = mate 	
MANDATORY_PATH = /usr/share/gconf/mate.mandatory.path 	
MATE_DESKTOP_SESSION_ID = this-is-deprecated 	
LESSOPEN = | /usr/bin/lesspipe %s 	
MDM_LANG = fr_BE.UTF-8 	
LOGNAME = tais 	
USER = tais 	
HOME = /home/tais 	
XDG_VTNR = 8 	
PATH = /usr/local/bin:/home/tais/BIN:/home/tais/bin:/home/tais/.local/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/home/tais/SRC/GRASS/grass_trunk/dist.x86_64-pc-linux-gnu/bin:/home/tais/SRC/GRASS/grass_trunk/dist.x86_64-pc-linux-gnu/script:/home/tais/SRC/GRASS/grass_trunk/dist.x86_64-pc-linux-gnu/lib:/home/tais/.grass7/addons:/home/tais/.grass7/addons/bin:/home/tais/.grass7/addons:/home/tais/.grass7/addons/scripts 	
CLICOLOR = 1 	
DISPLAY = :0.0 	
SSH_AGENT_PID = 2066 	
LANG = fr_BE.UTF-8 	
TERM = xterm-color 	
SHELL = /bin/bash 	
GIS_LOCK = $$ 	
XAUTHORITY = /home/tais/.Xauthority 	
SESSION_MANAGER = local/tais-HP-Z620-Workstation:@/tmp/.ICE-unix/1995,unix/tais-HP-Z620-Workstation:/tmp/.ICE[...]<jupyter_text>**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**   User inputs<jupyter_code>## Define a empty dictionnary for saving user inputs
user={}<jupyter_output><empty_output><jupyter_text>Here after:- Enter the path to the directory you want to use as "[GRASSDATA](https://grass.osgeo.org/programming7/loc_struct.png)". - Enter the name of the location in which you want to work and its projection information in [EPSG code](http://spatialreference.org/ref/epsg/) format. Please note that the GRASSDATA folder and locations will be automatically created if not existing yet. If the location name already exists, the projection information will not be used.  - Enter the name you want for the mapsets which will be used later for Unsupervised Segmentation Parameter Optimization (USPO), Segmentation and Classification steps.<jupyter_code>## Enter the path to GRASSDATA folder
user["gisdb"] = "/media/tais/My_Book_1/MAUPP/Traitement/Ouagadougou/Segmentation_fullAOI_localapproach/GRASSDATA"

## Enter the name of the location (existing or for a new one)
user["location"] = "Ouaga_32630"

## Enter the EPSG code for this location 
user["locationepsg"] = "32630"

## Enter the name of the mapset to use for segmentation
user["segmentation_mapsetname"] = "LOCAL_SEGMENT"

## Enter the name of the mapset to use for classification
user["classificationA_mapsetname"] = "CLASSIF"<jupyter_output><empty_output><jupyter_text>**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**   Define GRASSDATA folder and create GRASS' location and mapsets    Here after, the python script will check if GRASSDATA folder, locations and mapsets already exist. If not, they will be automatically created. **Import GRASS Python packages**<jupyter_code>## Import libraries needed to launch GRASS GIS in the jupyter notebook
import grass.script.setup as gsetup

## Import libraries needed to call GRASS using Python
import grass.script as grass<jupyter_output><empty_output><jupyter_text>**Define GRASSDATA folder and create location and mapsets**<jupyter_code>## Automatic creation of GRASSDATA folder
if os.path.exists(user["gisdb"]):
    print "GRASSDATA folder already exist" 
else: 
    os.makedirs(user["gisdb"]) 
    print "GRASSDATA folder created in "+user["gisdb"]
## Automatic creation of GRASS location is doesn't exist
if os.path.exists(os.path.join(user["gisdb"],user["location"])):
    print "Location "+user["location"]+" already exist" 
else : 
    grass.core.create_location(user["gisdb"], user["location"], epsg=user["locationepsg"], overwrite=False)
    print "Location "+user["location"]+" created"
### Automatic creation of GRASS GIS mapsets

## Import library for file copying 
import shutil

mapsetname=user["classificationA_mapsetname"]
if os.path.exists(os.path.join(user["gisdb"],user["location"],mapsetname)):
    if not os.path.exists(os.path.join(user["gisdb"],user["location"],mapsetname,'WIND')):
        print "WARNING: '"+mapsetname+"' mapset already exist, but a 'WIND' file is missing. Please solve this issue."
    else: print "'"+mapsetname+"' mapset already exist" 
else: 
    os.makedirs(os.path.join(user["gisdb"],user["location"],mapsetname))
    shutil.copy(os.path.join(user["gisdb"],user["location"],'PERMANENT','WIND'),os.path.join(user["gisdb"],user["location"],mapsetname,'WIND'))
    print "'"+mapsetname+"' mapset created in location '"+user["location"]+"'"
### Automatic creation of GRASS GIS mapsets

## Import library for file copying 
import shutil

mapsetname=user["segmentation_mapsetname"]
if os.path.exists(os.path.join(user["gisdb"],user["location"],mapsetname)):
    if not os.path.exists(os.path.join(user["gisdb"],user["location"],mapsetname,'WIND')):
        print "WARNING: '"+mapsetname+"' mapset already exist, but a 'WIND' file is missing. Please solve this issue."
    else: print "'"+mapsetname+"' mapset already exist" 
else: 
    os.makedirs(os.path.join(user["gisdb"],user["location"],mapsetname))
    shutil.copy(os.path.join(user["gisdb"],user["location"],'PERMANENT','WIND'),os.path.join(user["gisdb"],user["location"],mapsetname,'WIND'))
    print "'"+mapsetname+"' mapset created in location '"+user["location"]+"'"<jupyter_output>'LOCAL_SEGMENT' mapset already exist<jupyter_text>**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**   Define functions    This section of the notebook is dedicated to defining functions which will then be called later in the script. If you want to create your own functions, define them here.  Function for computing processing time The "print_processing_time" is used to calculate and display the processing time for various stages of the processing chain. At the beginning of each major step, the current time is stored in a new variable, using [time.time() function](https://docs.python.org/2/library/time.html). At the end of the stage in question, the "print_processing_time" function is called and takes as argument the name of this new variable containing the recorded time at the beginning of the stage, and an output message.<jupyter_code>## Import library for managing time in python
import time  

## Function "print_processing_time()" compute processing time and printing it.
# The argument "begintime" wait for a variable containing the begintime (result of time.time()) of the process for which to compute processing time.
# The argument "printmessage" wait for a string format with information about the process. 
def print_processing_time(begintime, printmessage):    
    endtime=time.time()           
    processtime=endtime-begintime
    remainingtime=processtime

    days=int((remainingtime)/86400)
    remainingtime-=(days*86400)
    hours=int((remainingtime)/3600)
    remainingtime-=(hours*3600)
    minutes=int((remainingtime)/60)
    remainingtime-=(minutes*60)
    seconds=round((remainingtime)%60,1)

    if processtime<60:
        finalprintmessage=str(printmessage)+str(seconds)+" seconds"
    elif processtime<3600:
        finalprintmessage=str(printmessage)+str(minutes)+" minutes and "+str(seconds)+" seconds"
    elif processtime<86400:
        finalprintmessage=str(printmessage)+str(hours)+" hours and "+str(minutes)+" minutes and "+str(seconds)+" seconds"
    elif processtime>=86400:
        finalprintmessage=str(printmessage)+str(days)+" days, "+str(hours)+" hours and "+str(minutes)+" minutes and "+str(seconds)+" seconds"
    
    return finalprintmessage<jupyter_output><empty_output><jupyter_text>Function for Postgres database vaccum<jupyter_code># Do a VACUUM on the current Postgresql database
def vacuum(db):
    old_isolation_level = db.isolation_level
    db.set_isolation_level(0)
    query = "VACUUM"
    cur.execute(query)
    db.set_isolation_level(old_isolation_level)<jupyter_output><empty_output><jupyter_text>**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**<jupyter_code>## Saving current time for processing time management
begintime_full=time.time()<jupyter_output><empty_output><jupyter_text>**-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-**  Set parameter for Postgresql database connection<jupyter_code># User for postgresql connexion
dbuser="tais"
# Password of user
dbpassword="<PASSWORD>"
# Host of database
host="localhost"
# Name of the new database
dbname="ouaga_fullaoi_localsegment"
# Set name of schema for objects statistics
stat_schema="statistics"
# Set name of schema for samples
sample_schema="samples"
# Set name of schema for classification optical only
classifA_schema="classif_A"
# Set name of schema for classification optical and nDSM
classifB_schema="classif_B"
# Set name of schema for classification SAR only
classifC_schema="classif_C"
# Set name of schema for classification SAR and optical
classifD_schema="classif_D"
# Set name of table with statistics of segmentobject_stats_table="object_stats_sar"s - FOR OPTICAL
object_stats_optical="object_stats_optical"
# Set name of table with statistics of segmentobject_stats_table="object_stats_sar"s - FOR SAR
object_stats_sar="object_stats_sar"
# Set name of table with all the samples
samples_labels="sample_labels"
# Set name of table with samples without outliers
samples_labels_ok="sample_labels_ok"
# Set name of table with extension of the AOI where SAR and optical exists
sar_aoi="sar_aoi"
# Set name of table with samples without outliers in SAR AOI
sample_sar="sample_sar"
# Set name of table with samples without outliers in SAR AOI - Test set
sample_test="sample_test"
# Set name of table with samples without outliers in SAR AOI - Training set
sample_training="sample_training"
# Set name of table with results of classification
classif="classif"
# Set name of table with results of classification and ground truth
groundtruth_classif="groundtruth_classif"<jupyter_output><empty_output><jupyter_text>Select samples in segmentation layer     **Launch GRASS GIS working session**<jupyter_code>### Automatic creation of GRASS GIS mapsets

## Import library for file copying 
import shutil

## Set the name of the mapset in which to work
mapsetname=classifC_schema

if os.path.exists(os.path.join(user["gisdb"],user["location"],mapsetname)):
    if not os.path.exists(os.path.join(user["gisdb"],user["location"],mapsetname,'WIND')):
        print "WARNING: '"+mapsetname+"' mapset already exist, but a 'WIND' file is missing. Please solve this issue."
    else: print "'"+mapsetname+"' mapset already exist" 
else: 
    os.makedirs(os.path.join(user["gisdb"],user["location"],mapsetname))
    shutil.copy(os.path.join(user["gisdb"],user["location"],'PERMANENT','WIND'),os.path.join(user["gisdb"],user["location"],mapsetname,'WIND'))
    print "'"+mapsetname+"' mapset created in location '"+user["location"]+"'"
## Launch GRASS GIS working session in the mapset
if os.path.exists(os.path.join(user["gisdb"],user["location"],mapsetname)):
    gsetup.init(os.environ['GISBASE'], user["gisdb"], user["location"], mapsetname)
    print "You are now working in mapset '"+mapsetname+"'" 
else: 
    print "'"+mapsetname+"' mapset doesn't exists in "+user["gisdb"]<jupyter_output>You are now working in mapset 'classif_C'<jupyter_text>Create new schema<jupyter_code>from psycopg2.extensions import ISOLATION_LEVEL_AUTOCOMMIT

# Connect to postgres database
db=None
db=pg.connect(dbname=dbname, user='tais', password='<PASSWORD>', host='localhost')

# Allow to create a new database
db.set_isolation_level(ISOLATION_LEVEL_AUTOCOMMIT)

# Execute the CREATE DATABASE query
cur=db.cursor()
#cur.execute('DROP SCHEMA IF EXISTS '+classifC_schema+' CASCADE') #Comment this to avoid deleting existing DB
try:
    cur.execute('CREATE SCHEMA '+classifC_schema)
except Exception as e:
    print ("Exception occured : "+str(e))
cur.close()
db.close()<jupyter_output>Exception occured : ERREUR:  le schéma « classif_c » existe déjà<jupyter_text>Define the folder where to save the results and create it if necessary In the next cell, please adapt the path to the directory where you want to save the .csv output of i.segment.uspo.<jupyter_code>general_path="/media/tais/My_Book_1/MAUPP/Traitement/Ouagadougou/Segmentation_fullAOI_localapproach/Results/CLASSIF"
resultfolder=os.path.join(general_path, classifC_schema)
## Create the folder if does not exists
if not os.path.exists(resultfolder):
    os.makedirs(resultfolder)
    print "Folder '"+resultfolder+"' created"
outputfolder=os.path.join(resultfolder, "classification")
## Create the folder if does not exists
if not os.path.exists(outputfolder):
    os.makedirs(outputfolder)
    print "Folder '"+outputfolder+"' created"<jupyter_output><empty_output><jupyter_text>Create a README.txt file informing about the classification methodology<jupyter_code>## Set-up a message
message=""
message+=classifC_schema+" : Classification using SAR only"
## Write it in the .txt file
f=open(os.path.join(outputfolder,"readme.txt"),'w')
f.write(message)
f.close()<jupyter_output><empty_output><jupyter_text>Display test and training samples<jupyter_code># Connect to an existing database
db=pg.connect(database=dbname, user=dbuser, password=<PASSWORD>, host=host)
# Open a cursor to perform database operations
cur=db.cursor()
# Query to find the number of row in the test samples
query="SELECT class, min(class_num) as class_num, count(*) FROM "+sample_schema+"."+sample_test+" \
GROUP BY class \
ORDER BY max(class_num) ASC"
# Execute query through panda
df_countperclass=pd.read_sql(query, db)
# Show dataframe
df_countperclass.head(5000)
# Print the total number in sample 
print "Test samples counts "+str(sum(df_countperclass['count']))+" points."
# Query to find the number of row in the sample table
query="SELECT class, min(class_num) as class_num, count(*) FROM "+sample_schema+"."+sample_training+" \
GROUP BY class \
ORDER BY max(class_num) ASC"
# Execute query through panda
df_countperclass=pd.read_sql(query, db)
# Show dataframe
df_countperclass.head(5000)
# Print the total number in sample 
print "Training samples counts "+str(sum(df_countperclass['count']))+" points."<jupyter_output>Training samples counts 1704 points.<jupyter_text>Create csv with statistics for training<jupyter_code># Drop table if exists:
cur.execute("DROP TABLE IF EXISTS "+classifC_schema+"."+sample_training)
# Make the changes to the database persistent
db.commit()
# Create table with 
query="CREATE TABLE "+classifC_schema+"."+sample_training+" AS(\
WITH shape_stats AS(\
SELECT a.seg_id, a.class_num, b.area, b.perimeter, b.compact_circle, b.compact_square, b.fd \
FROM "+sample_schema+"."+sample_training+" AS a \
INNER JOIN "+stat_schema+"."+object_stats_optical+" AS b \
ON a.seg_id=b.cat)\
\
SELECT a.*, c.* \
FROM shape_stats AS a \
INNER JOIN "+stat_schema+"."+object_stats_sar+" AS c \
ON a.seg_id=c.cat)"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
# Drop some non wanted columns
query="ALTER TABLE "+classifC_schema+"."+sample_training+" DROP COLUMN cat"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
# Query to find the number of row in the sample table
query="SELECT * FROM "+classifC_schema+"."+sample_training
# Execute query through panda
df_countperclass=pd.read_sql(query, db)
# Show dataframe
df_countperclass.head(10)<jupyter_output><empty_output><jupyter_text>Export training sample as .csv for archive<jupyter_code>outputfolder
## Define the path to the .csv output
sample_training_csv=os.path.join(outputfolder,"sample_training.csv")
# Connect to an existing database
db=pg.connect(database=dbname, user=dbuser, password=<PASSWORD>, host=host)
# Open a cursor to perform database operations
cur=db.cursor()
#### Export as .csv
# Query
query="COPY "+classifC_schema+"."+sample_training+" TO '"+sample_training_csv+"' DELIMITER ',' CSV HEADER"
# Execute the CREATE TABLE query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
# Close cursor and communication with the database
cur.close()
db.close()<jupyter_output><empty_output><jupyter_text>Feature selection      Import and load extension rpy2<jupyter_code>import rpy2
%load_ext rpy2.ipython<jupyter_output><empty_output><jupyter_text>Import training set in R<jupyter_code>%%R -i sample_training_csv -o training
# Import training
training <- read.csv(sample_training_csv, sep=",", header=TRUE, row.names=1)
training$class_num <- as.factor(training$class_num)
training.head()
%%R -o colname
colname<-colnames(training)
# Create a list with the column name
print list(colname)
# Set the name of the colum with class
class_column='class_num'
# Set the name of the colum with the first feature to be used
first_feature='area'
# Set the name of the colum with the last feature to be used
last_feature='w11_sm_perc_90'
%%R -i class_column,first_feature,last_feature -o features

# Save the index of the column
classnumindex=match(c(class_column),colname)
# Save the index where object's features to be used start
startindex=match(c(first_feature),colname)
# Save the index where object's features to be used stop
stopindex=match(c(last_feature),colname)

# Create dataframe with only columns to be used features
temp1<-as.data.frame(training[classnumindex])
temp2<-as.data.frame(training[startindex:stopindex])
merged_df=merge(temp1,temp2, by=0, all=TRUE)
features<-transform(merged_df, row.names=Row.names, Row.names=NULL)
# Set class column as factor
features$class_num <- as.factor(features$class_num)
features.head()
%%R 
# Instal library
if(!is.element('parallel', installed.packages()[,1])){
    cat('\n\nInstalling parallel package from CRAN')
    chooseCRANmirror(ind=1)
    install.packages('parallel')}

if(!is.element('doParallel', installed.packages()[,1])){
    cat('\n\nInstalling doParallel package from CRAN')
    chooseCRANmirror(ind=1)
    install.packages('doParallel')}

if(!is.element('VSURF', installed.packages()[,1])){
    cat('\n\nInstalling VSURF package from CRAN')
    chooseCRANmirror(ind=1)
    install.packages('VSURF')}

library(parallel)
require(parallel)
library(doParallel)
require(doParallel)
library(VSURF)
require(VSURF)
%%R -i class_column,first_feature,last_feature -o Surf

# Save column names
colname<-colnames(features)
# Save the index of the column
classnumindex<-match(c(class_column),colname)
# Save the index where object's features to be used start
startindex<-match(c(first_feature),colname)
# Save the index where object's features to be used stop
stopindex<-match(c(last_feature),colname)

# Set number of cores to use
usedcores=detectCores()-5
registerDoParallel(usedcores)
# Feature Selection using VSURF
Surf=VSURF(features[startindex:stopindex],features$class_num, parallel=TRUE, ncores=usedcores)
## Set the path to the output .txt file with summary of VSURF
output_summary=os.path.join(outputfolder,"VSURF_summary.txt")
%%R -i output_summary

# Print and save summary of VSRUF
summaryVSRUF<-summary(Surf)
output<-capture.output(summary(Surf))
print (summaryVSRUF)
cat(output,file=output_summary,sep="\n")
## Set the path to the output csv with results of VSURF
output_interp=os.path.join(outputfolder,"VSURF_interp.csv")
output_pred=os.path.join(outputfolder,"VSURF_pred.csv")
%%R -i output_interp,output_pred -o Surf_interp_features,Surf_pred_features

# Save interpretation step and prediction step results
Surf_interp<-Surf$varselect.interp
Surf_pred<-Surf$varselect.pred

# Save name of feature in a list ranked with the same index as VSURF output
colname<-colnames(features[startindex:stopindex])

# Declare empty variables of type 'character'
Surf_interp_features <- character()

# Loop on indexes of features selecteds by VSURF to find the corresponding feature name
count<-0
for (x in Surf_interp){
  count<-count+1
  Surf_interp_features[count]<-colname[x]
}

# Write CSV
write.csv(Surf_interp_features, file=output_interp, row.names=TRUE)

# Declare empty variables of type 'character'
Surf_pred_features <- character()

# Loop on indexes of features selecteds by VSURF to find the corresponding feature name
count<-0
for (x in Surf_pred){
  count<-count+1
  Surf_pred_features[count]<-colname[x]
}

# Write CSV
write.csv(Surf_pred_features, file=output_pred, row.names=TRUE)
## Print features selected by VSURF at interpretation step
print "Interpretation results"
for i in range(len(Surf_interp_features)):
    print str(i+1)+"    "+Surf_interp_features[i] 
## Print features selected by VSURF at interpretation step
print "Prediction results"
for i in range(len(Surf_pred_features)):
    print str(i+1)+"    "+Surf_pred_features[i] 
## Set the path to the output plot (PDF)
outputgraph=os.path.join(outputfolder,"VSURF_opt_plot.pdf")
%%R -i outputgraph

# Export the plot in pdf
VSURF_plot_export=outputgraph
pdf(VSURF_plot_export)
plot(Surf, var.names = FALSE,
     nvar.interp = length(Surf_pred))
dev.off()
## Display graphs
import wand
from wand.image import Image as WImage
img = WImage(filename=outputgraph)
img<jupyter_output><empty_output><jupyter_text>Classification using Random Forest      Prepare .csv for classification (training set and object to classify)<jupyter_code># Connect to an existing database
db=pg.connect(database=dbname, user=dbuser, password=<PASSWORD>, host=host)
# Open a cursor to perform database operations
cur=db.cursor()<jupyter_output><empty_output><jupyter_text>Export training set as .csv<jupyter_code>### Set if all feature have to be used (vsurfonly=False) 
### or only those resulting from prediction step of VSURF (vsurfonly=True)
vsurfonly=True
## Save list of features to be used 
if vsurfonly:
    features_for_classif=list(Surf_pred_features)
else:
    features_for_classif=list(features.columns.values)[1:]
    
## Split features in to a list for Geometrical and Optical features and another for SAR features
features_for_classif_optical=[]
features_for_classif_sar=[]
for f in features_for_classif:
    if f.split("_")[0] in ("w3","w5","w7","w11","w13","w15","w17"):
        features_for_classif_sar.append(f)
    else:
        features_for_classif_optical.append(f)
print "Optical and geometric features to be used :\n"+','.join(features_for_classif_optical)
print "\n"
print "SAR features to be used :\n"+','.join(features_for_classif_sar)
print "\n"
print "All features to be used :\n"+','.join(features_for_classif)
## Define the path to the .csv
training_csv=os.path.join(outputfolder,"training_csv.csv")
#### Export as .csv
# Query
query="COPY (\
SELECT seg_id, class_num, "+", ".join(features_for_classif)+" \
FROM "+classifC_schema+"."+sample_training+") TO '"+training_csv+"' DELIMITER ',' CSV HEADER"
# Execute the CREATE TABLE query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()<jupyter_output><empty_output><jupyter_text>Evaluate number of objects to be classified and number of loops to create<jupyter_code>## Evalutate the number of row in 'object_stat' table
query="SELECT reltuples::bigint AS estimated_object \
FROM pg_class \
WHERE oid='"+stat_schema+"."+object_stats_optical+"'::regclass"  
# Execute query through panda
df=pd.read_sql(query, db)
# Save estimated number of objects
estimate=list(df['estimated_object'])[0]
# Print
print "The table contains "+str(estimate)+" rows (estimate)."
# Define number of object to classify on each loop
nbobjloop=1000000
# Add 10% of estimated number of rows to be sure to well cover the full dataset 
remaining_rows=int(estimate*1.05)
# Define number of loops for classification
import math
loops=int(math.ceil(remaining_rows/(nbobjloop*1.0)))
## Print what is going to happend
print "Classification will be made on "+str(loops)+" loop(s)"<jupyter_output>Classification will be made on 17 loop(s)<jupyter_text>Create mutliple .csv with the statistics of segments  Creating a new temporary table containing objects to be classified with all the statistics needed for classification<jupyter_code># Query to create the temporary table
query="CREATE TABLE "+classifC_schema+".temp1 AS(\
WITH \
optistats AS(\
SELECT cat, "+", ".join(features_for_classif_optical)+" \
FROM "+stat_schema+"."+object_stats_optical+"), \
\
sarstats AS(\
SELECT cat, "+", ".join(features_for_classif_sar)+" \
FROM "+stat_schema+"."+object_stats_sar+"),\
\
allstats AS(\
SELECT a.cat AS acat, a.*, b.* \
FROM optistats AS a \
INNER JOIN sarstats AS b \
ON a.cat=b.cat)\
\
SELECT acat AS cat, "+", ".join(features_for_classif)+" \
FROM allstats ORDER BY cat)"
# Execute the query 
cur.execute(query)
db.commit()
# Define starting offset and limit (for the firsth loop)
offset=0
# Declare an empty list containing the paths to .csv
list_objstat_csv=[]
# Export object statistics on multiple .csv files
for loop in range(loops):
    ## Define the path to the .csv
    objstat_csv=os.path.join(outputfolder,"objects_stats_"+str(loop+1)+".csv")
    list_objstat_csv.append(objstat_csv)        
    # Query to export as .csv
    query="COPY (SELECT * FROM "+classifC_schema+".temp1 ORDER BY cat \
    OFFSET "+str(offset)+" LIMIT "+str(nbobjloop)+") TO '"+objstat_csv+"' DELIMITER ',' CSV HEADER"
    # Execute the query 
    cur.execute(query)
    db.commit()
    # Print
    print "Rows between "+str(offset)+" and "+str(offset+nbobjloop)+" exported in .csv file '"+objstat_csv+"'"
    # Update offset and limit 
    offset+=nbobjloop
# Drop table if exists:
cur.execute("DROP TABLE IF EXISTS "+classifC_schema+".temp1")
# Make the changes to the database persistent
db.commit()
# Close cursor and communication with the database
cur.close()
db.close()
%%R 

# Install package
if(!is.element('caret', installed.packages()[,1])){
  cat('\n\nInstalling caret package from CRAN')
  chooseCRANmirror(ind=1)
  install.packages('caret')}
# Install package
if(!is.element('randomForest', installed.packages()[,1])){
  cat('\n\nInstalling randomForest package from CRAN')
  chooseCRANmirror(ind=1)
  install.packages('randomForest')}
# Install package
if(!is.element('pROC', installed.packages()[,1])){
  cat('\n\nInstalling pROC package from CRAN')
  chooseCRANmirror(ind=1)
  install.packages('pROC')}
# Install package
if(!is.element('e1071', installed.packages()[,1])){
  cat('\n\nInstalling e1071 package from CRAN')
  chooseCRANmirror(ind=1)
  install.packages('e1071')}
%%R 

# Load libraries
library(parallel)
require(parallel)
library(doParallel)
require(doParallel)
library(caret)
require(caret)
library(randomForest)
require(randomForest)
library(e1071)
require(e1071)<jupyter_output>/usr/local/lib/python2.7/dist-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: Le chargement a nécessité le package : lattice

  warnings.warn(x, RRuntimeWarning)
/usr/local/lib/python2.7/dist-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: Le chargement a nécessité le package : ggplot2

  warnings.warn(x, RRuntimeWarning)
/usr/local/lib/python2.7/dist-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: randomForest 4.6-12

  warnings.warn(x, RRuntimeWarning)
/usr/local/lib/python2.7/dist-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: Type rfNews() to see new features/changes/bug fixes.

  warnings.warn(x, RRuntimeWarning)
/usr/local/lib/python2.7/dist-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: 
Attachement du package : ‘randomForest’


  warnings.warn(x, RRuntimeWarning)
/usr/local/lib/python2.7/dist-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: The following object is masked from ‘package:ggplot2’:

    margin


  w[...]<jupyter_text>Train Random Forest model<jupyter_code>## Set the path to the output of random forest classification
output_accuracy=os.path.join(outputfolder,"accuracy_RF.csv")
%%R -i training_csv,output_accuracy

# Set number of cores to use
usedcores=detectCores()-5
registerDoParallel(usedcores)

# Import training
training <- read.csv(training_csv, sep=",", header=TRUE, row.names=1)
# Define factor
training$class_num <- as.factor(training$class_num)

# Cross-validation setting
MyFolds.cv <- createMultiFolds(training$class_num, k=5, times=10)
MyControl.cv <- trainControl(method='repeatedCV', index=MyFolds.cv, allowParallel = TRUE)

# Train Random Forest
rfModel <- train(class_num~.,training,method='rf', trControl=MyControl.cv,tuneLength=10)
resamps.cv <- rfModel$resample
accuracy_means <- mean(resamps.cv$Accuracy)
kappa_means <- mean(resamps.cv$Kappa)
df_means <- data.frame(method='rf',accuracy=accuracy_means, kappa=kappa_means)
write.csv(df_means, output_accuracy, row.names=FALSE, quote=FALSE)
#### Show mean accuracy results from cross-validation for tuning
## Import .csv file
accuracy=pd.read_csv(output_accuracy, sep=',',header=0)
## Display table
accuracy.head(15)
## Set the path to the output of random forest classification
output_varimp=os.path.join(outputfolder,"VariablesImportance_RF.pdf")
output_rfmodel=os.path.join(outputfolder,"RF_model.txt")
output_cv=os.path.join(outputfolder,"RF_cv.pdf")
%%R -i output_varimp,output_rfmodel,output_cv

# Plot variable importance
library(pROC)
importance <- varImp(rfModel, scale=FALSE)
pdf(output_varimp, width = 11, height = ncol(features)/6 )
print(plot(importance))
dev.off()

# Show final model
tmp<-rfModel$finalModel
sink(output_rfmodel)
print(tmp)
sink()

# Plot cross validation tuning results
pdf(output_cv)
print(plot(rfModel))
dev.off()
#### Import classifiers tuning parameters and confusion matrix
## Open file
classifier_runs = open(output_rfmodel, 'r')  
## Read file
print classifier_runs.read()<jupyter_output>Call:
 randomForest(x = x, y = y, mtry = param$mtry) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 2

        OOB estimate of  error rate: 61.56%
Confusion matrix:
    11 13 14 21 22 31 32 33 34 41 class.error
11 283 31  0  5  0 10  0  0 26  7   0.2182320
13  69 50  0  0  1  4  0  0  3  0   0.6062992
14   0  1 28 18  8  0 15  3 21  5   0.7171717
21  15  0 16 81  4  1  7  7 20  5   0.4807692
22  10  2  3 11 26  8 20 12 24  6   0.7868852
31  71 16  4  2  5 17  5  2 22  2   0.8835616
32  12  0  7 11 13  5 40 15 41  9   0.7385621
33   5  0  6 10  8  3 25 18 49  7   0.8625954
34  54  8 12 13  9 10 31 19 85 19   0.6730769
41  33  2  4  5  5  9 11 13 39 27   0.8175676<jupyter_text>The following cell display pdf figures. Please first install the **Wand** package using the folowwing command in the terminal : ``` sudo pip instal Wand ```<jupyter_code>## Display graphs
import wand
from wand.image import Image as WImage
img = WImage(filename=output_varimp)
img
## Display graphs
import wand
from wand.image import Image as WImage
img = WImage(filename=output_cv)
img<jupyter_output><empty_output><jupyter_text>Predict classes for objects<jupyter_code>%%R -i list_objstat_csv

# Classification using loop
for(i in 1:length(list_objstat_csv)) {
    filepath<-list_objstat_csv[i]
    # Import object statistics as dataframe
    features <- read.csv(filepath, sep=",", header=TRUE, row.names=1)
    # Predict class 
    predicted <- data.frame(predict(rfModel, features))
    name<-paste('resultsdf',i,sep='')
    assign(name,data.frame(id=rownames(features), predicted))
    rm(features)
    rm(predicted)
    gc()
}
%%R

# Merge predictions in one single dataframe
listdf<-lapply(ls(pattern = "resultsdf*"), get)
rf_predictions<-do.call(rbind,listdf)<jupyter_output><empty_output><jupyter_text>Save training set and predictions in .csv files<jupyter_code>## Define outputfile for .csv with training data for R
trainingset_R=os.path.join(outputfolder,"RF_trainingset_R.csv")
## Define outputfile for .csv with predictions of classification
classif_results=os.path.join(outputfolder,"RF_classif_results_R.csv")
%%R -i trainingset_R,classif_results
write.csv(training, file=trainingset_R, row.names=FALSE, quote=FALSE)
write.csv(rf_predictions, file=classif_results, row.names=FALSE, quote=FALSE)
%%R 

# Remove variables not needed anymore
rm(training)
gc()

rm(list=ls(pattern = "resultsdf*"))
rm(listdf)
gc()<jupyter_output><empty_output><jupyter_text>Delete csv files with object statistics for each loop of classification<jupyter_code>for csv in list_objstat_csv:
    os.remove(csv)<jupyter_output><empty_output><jupyter_text>Import classification results in postgis<jupyter_code># Connect to an existing database
db=pg.connect(database=dbname, user=dbuser, password=<PASSWORD>, host=host)
# Open a cursor to perform database operations
cur=db.cursor()<jupyter_output><empty_output><jupyter_text>Import predictions of Random forest at level 2<jupyter_code># Drop table if exists:
cur.execute("DROP TABLE IF EXISTS "+classifC_schema+"."+classif)
# Make the changes to the database persistent
db.commit()

# Create new table
query="CREATE TABLE "+classifC_schema+"."+classif+" (seg_id integer PRIMARY KEY, rf_pred_l2 integer)"
# Execute the CREATE TABLE query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()

# Create query for copy data from csv
query="COPY "+classifC_schema+"."+classif+" FROM '"+str(classif_results)+"' HEADER DELIMITER ',' CSV" 
# Execute the COPY FROM CSV query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()<jupyter_output><empty_output><jupyter_text>Merge classes of level 2 to reach level 2b<jupyter_code>## Add column for prediction at level 2b
query="ALTER TABLE "+classifC_schema+"."+classif+" ADD COLUMN rf_pred_l2b integer"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
# Update the value
query="UPDATE "+classifC_schema+"."+classif+" SET \
rf_pred_l2b = '20' WHERE  rf_pred_l2 IN ('21','22')"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
# Update the value
query="UPDATE "+classifC_schema+"."+classif+" SET \
rf_pred_l2b = '30' WHERE  rf_pred_l2 IN ('32','33','34')"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
# Update the value
query="UPDATE "+classifC_schema+"."+classif+" SET \
rf_pred_l2b = rf_pred_l2 WHERE rf_pred_l2b IS NULL"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()<jupyter_output><empty_output><jupyter_text>Merge classes of level 2 to reach level 1<jupyter_code>## Add column for prediction at level 1
query="ALTER TABLE "+classifC_schema+"."+classif+" ADD COLUMN rf_pred_l1 integer"
# Execute the query 
cur.execute(query)
## Update column for prediction and ground truth at level 2 
query="UPDATE "+classifC_schema+"."+classif+" SET rf_pred_l1=rf_pred_l2/10"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()<jupyter_output><empty_output><jupyter_text>Export to csv for archive<jupyter_code>## Define path of .txt file for saving prediction at different level
outputfile=classif_results.split(".")[0]+'_diff_levels.csv'
#### Export as .csv
# Query
query="COPY "+classifC_schema+"."+classif+" TO '"+outputfile+"' DELIMITER ',' CSV HEADER"
# Execute the CREATE TABLE query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
# Close cursor and communication with the database
cur.close()
db.close()<jupyter_output><empty_output><jupyter_text>Accuracy assesment      Change the outputfolder<jupyter_code>## Folder in which save processing time output
outputfolder=os.path.join(resultfolder,"accuracy_assess")

## Create the folder if does not exists
if not os.path.exists(outputfolder):
    os.makedirs(outputfolder)
    print "Folder '"+outputfolder+"' created"<jupyter_output><empty_output><jupyter_text>Join sample test and classifier predictions in postgresql<jupyter_code># Connect to an existing database
db=pg.connect(database=dbname, user=dbuser, password=<PASSWORD>, host=host)
# Open a cursor to perform database operations
cur=db.cursor()
# Build a query to drop view if exists
query="DROP TABLE IF EXISTS "+classifC_schema+"."+groundtruth_classif
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
#### Join sample_test and classification results
query="CREATE TABLE "+classifC_schema+"."+groundtruth_classif+" AS (\
SELECT a.cat_point, a.id, \
a.class AS level2_label, \
a.class_b AS level2b_label, \
a.class_l1 AS level1_label, \
a.class_num::integer AS level2_groundtr, \
a.class_num_b::integer AS level2b_groundtr, \
a.class_num_l1::integer AS level1_groundtr, \
c.rf_pred_l2 AS rf_level2, \
c.rf_pred_l2b AS rf_level2b, \
c.rf_pred_l1 AS rf_level1 \
FROM "+sample_schema+"."+sample_test+" AS a \
LEFT JOIN "+classifC_schema+"."+classif+" AS c ON a.seg_id = c.seg_id)"
# Execute the query 
cur.execute(query)
# Make the changes to the database persistent
db.commit()
#### Save ground truth, optical_predic, opticalsar_predic as lists, for Level-1 and Level-2
# Query 
query="SELECT * FROM "+classifC_schema+"."+groundtruth_classif+" ORDER BY id"
# Execute query through panda
df=pd.read_sql(query, db)

# Show dataframe
df.head(20)
# Save ground truth - Level2
groundtruth_L2=list(df['level2_groundtr'])
# Save ground truth - Level2b
groundtruth_L2b=list(df['level2b_groundtr'])
# Save ground truth - Level1
groundtruth_L1=list(df['level1_groundtr'])

# Save predictions - level2
prediction_L2=list(df['rf_level2'])
# Save predictions - level2b
prediction_L2b=list(df['rf_level2b'])
# Save predictions - level1
prediction_L1=list(df['rf_level1'])
# Display number of objects in test sample
nrows=len(groundtruth_L2)
print nrows
## Zip the list together - Level 2
ziped_results=zip(groundtruth_L2,prediction_L2)
## Sort
ziped_results.sort()
## Unzip
groundtruth_L2,prediction_L2=zip(*ziped_results)
## Zip the list together - Level 2b
ziped_results=zip(groundtruth_L2b,prediction_L2b)
## Sort
ziped_results.sort()
## Unzip
groundtruth_L2b,prediction_L2b=zip(*ziped_results)
## Zip the list together - Level 1
ziped_results=zip(groundtruth_L1,prediction_L1)
## Sort
ziped_results.sort()
## Unzip
groundtruth_L1,prediction_L1=zip(*ziped_results)
## Check if lists contain th same distinct values
tmp1=list(set(groundtruth_L2))
tmp1.sort()
tmp2=list(set(prediction_L2))
tmp2.sort()
if tmp1 != tmp2:
    sys.exit('WARNING: Lists contain different distinct values. Please check before continue')
    
## Check if lists contain th same distinct values
tmp1=list(set(groundtruth_L2b))
tmp1.sort()
tmp2=list(set(prediction_L2b))
tmp2.sort()
if tmp1 != tmp2:
    sys.exit('WARNING: Lists contain different distinct values. Please check before continue')

## Check if lists contain th same distinct values
tmp1=list(set(groundtruth_L1))
tmp1.sort()
tmp2=list(set(prediction_L1))
tmp2.sort()
if tmp1 != tmp2:
    sys.exit('WARNING: Lists contain different distinct values. Please check before continue')
#### Save ground truth, optical_predic, opticalsar_predic as lists, for Level-1 and Level-2
# Query 
query="SELECT DISTINCT level2_groundtr, level2_label \
FROM "+classifC_schema+"."+groundtruth_classif+" ORDER BY level2_groundtr"
# Execute query through panda
df=pd.read_sql(query, db)
# Save ground truth classes labels - Level2
classes_L2=list(df['level2_label'])

# Query 
query="SELECT DISTINCT level2b_groundtr, level2b_label \
FROM "+classifC_schema+"."+groundtruth_classif+" ORDER BY level2b_groundtr"
# Execute query through panda
df=pd.read_sql(query, db)
# Save ground truth classes labels - Level2
classes_L2b=list(df['level2b_label'])

# Query 
query="SELECT DISTINCT level1_groundtr, level1_label \
FROM "+classifC_schema+"."+groundtruth_classif+" ORDER BY level1_groundtr"
# Execute query through panda
df=pd.read_sql(query, db)
# Save ground truth classes labels - Level2
classes_L1=list(df['level1_label'])
print "Level 2 classes:\n\n"+'\n'.join(classes_L2)
print "Level 2 classes:\n\n"+'\n'.join(classes_L2b)
print "Level 2 classes:\n\n"+'\n'.join(classes_L1)<jupyter_output>Level 2 classes:

Artificial surfaces
Natural material surfaces
Vegetation
Water<jupyter_text>Confusion matrix   The classification performance evaluation is made mainly using [Scikit-learn](http://scikit-learn.org/stable/modules/model_evaluation.html).The trick to plot the confusion matrix was found [here](http://scikit-learn.org/stable/auto_examples/model_selection/plot_confusion_matrix.htmlsphx-glr-auto-examples-model-selection-plot-confusion-matrix-py).<jupyter_code>## Import libraries
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix
import itertools
import numpy as np
%matplotlib inline  
def plot_confusion_matrix(cm, classes,
                          normalize=False,
                          title='Confusion matrix',
                          cmap=plt.cm.Blues):
    """
    This function prints and plots the confusion matrix.
    Normalization can be applied by setting `normalize=True`.
    """
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    tick_marks = np.arange(len(classes))
    plt.xticks(tick_marks, classes, ha='right', rotation=45)
    plt.yticks(tick_marks, classes)

    if normalize:
        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]
        print("Normalized confusion matrix")
    else:
        print('Confusion matrix, without normalization')

    #print(cm)

    thresh = cm.max() / 2.
    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
        if normalize:
            plt.text(j, i, round(cm[i, j],2),
            horizontalalignment="center",
            color="white" if cm[i, j] > thresh else "black")
        else:
            plt.text(j, i, cm[i, j],
            horizontalalignment="center",
            color="white" if cm[i, j] > thresh else "black")

    plt.tight_layout()
    plt.ylabel('True label')
    plt.xlabel('Predicted label')<jupyter_output><empty_output><jupyter_text>Level 2b<jupyter_code># Compute confusion matrix
cnf_matrix_L2=confusion_matrix(groundtruth_L2, prediction_L2)
## Set the path to the output
output_rowconfmat=os.path.join(outputfolder,"rowconfusionmatrix_L2.txt")

## Export the row confusion matrix
numpy.savetxt(output_rowconfmat, cnf_matrix_L2.astype(np.int), fmt='%d', delimiter=",")
# Plot non-normalized confusion matrix
fig_cm=plt.figure(figsize=(15,10))
plot_confusion_matrix(cnf_matrix_L2, classes=classes_L2,
                      title='Confusion matrix, without normalization')
# Plot normalized confusion matrix
fig_cm_normal=plt.figure(figsize=(15,10))
plot_confusion_matrix(cnf_matrix_L2, classes=classes_L2, normalize=True,
                      title='Normalized confusion matrix')

plt.show()
## Set the path to the output
output_confmat_pdf=os.path.join(outputfolder,"confusionmatrix_L2.pdf")
output_confmatA_png=os.path.join(outputfolder,"confusionmatrixA_L2.png")
output_confmatB_png=os.path.join(outputfolder,"confusionmatrixB_L2.png")
# Export in PDF
from matplotlib.backends.backend_pdf import PdfPages
pp=PdfPages(output_confmat_pdf)
pp.savefig(fig_cm)
pp.savefig(fig_cm_normal)
pp.close()

# Export in PNG
fig_cm.savefig(output_confmatA_png, format='png', dpi=300)
fig_cm_normal.savefig(output_confmatB_png, format='png', dpi=300)<jupyter_output><empty_output><jupyter_text>Level 2<jupyter_code># Compute confusion matrix
cnf_matrix_L2b=confusion_matrix(groundtruth_L2b, prediction_L2b)
## Set the path to the output
output_rowconfmat=os.path.join(outputfolder,"rowconfusionmatrix_L2b.txt")

## Export the row confusion matrix
numpy.savetxt(output_rowconfmat, cnf_matrix_L2b.astype(np.int), fmt='%d', delimiter=",")
# Plot non-normalized confusion matrix
fig_cm=plt.figure(figsize=(15,10))
plot_confusion_matrix(cnf_matrix_L2b, classes=classes_L2b,
                      title='Confusion matrix, without normalization')
# Plot normalized confusion matrix
fig_cm_normal=plt.figure(figsize=(15,10))
plot_confusion_matrix(cnf_matrix_L2b, classes=classes_L2b, normalize=True,
                      title='Normalized confusion matrix')

plt.show()
## Set the path to the output
output_confmat_pdf=os.path.join(outputfolder,"confusionmatrix_L2b.pdf")
output_confmatA_png=os.path.join(outputfolder,"confusionmatrixA_L2b.png")
output_confmatB_png=os.path.join(outputfolder,"confusionmatrixB_L2b.png")
# Export in PDF
from matplotlib.backends.backend_pdf import PdfPages
pp=PdfPages(output_confmat_pdf)
pp.savefig(fig_cm)
pp.savefig(fig_cm_normal)
pp.close()

# Export in PNG
fig_cm.savefig(output_confmatA_png, format='png', dpi=300)
fig_cm_normal.savefig(output_confmatB_png, format='png', dpi=300)<jupyter_output><empty_output><jupyter_text>Level 1<jupyter_code># Compute confusion matrix
cnf_matrix_L1=confusion_matrix(groundtruth_L1, prediction_L1)
## Set the path to the output
output_rowconfmat=os.path.join(outputfolder,"rowconfusionmatrix_L1.txt")

## Export the row confusion matrix
numpy.savetxt(output_rowconfmat, cnf_matrix_L1.astype(np.int), fmt='%d', delimiter=",")
# Plot non-normalized confusion matrix
fig_cm=plt.figure(figsize=(10,7))
plot_confusion_matrix(cnf_matrix_L1, classes=classes_L1,
                      title='Confusion matrix, without normalization')
# Plot normalized confusion matrix
fig_cm_normal=plt.figure(figsize=(10,7))
plot_confusion_matrix(cnf_matrix_L1, classes=classes_L1, normalize=True,
                      title='Normalized confusion matrix')

plt.show()
## Set the path to the output
output_confmat_pdf=os.path.join(outputfolder,"confusionmatrix_L1.pdf")
output_confmatA_png=os.path.join(outputfolder,"confusionmatrixA_L1.png")
output_confmatB_png=os.path.join(outputfolder,"confusionmatrixB_L1.png")
# Export in PDF
from matplotlib.backends.backend_pdf import PdfPages
pp=PdfPages(output_confmat_pdf)
pp.savefig(fig_cm)
pp.savefig(fig_cm_normal)
pp.close()

# Export in PNG
fig_cm.savefig(output_confmatA_png, format='png', dpi=300)
fig_cm_normal.savefig(output_confmatB_png, format='png', dpi=300)<jupyter_output><empty_output><jupyter_text>Accuracy measures    Classification report TODO: Find a way to compute comission and omission errore / User's Producers' accuracies<jupyter_code>## Import libraries
from sklearn.metrics import accuracy_score
from sklearn.metrics import cohen_kappa_score
from sklearn.metrics import classification_report

from sklearn.metrics import f1_score<jupyter_output><empty_output><jupyter_text>Level 2<jupyter_code>## Set the path to the output
output=os.path.join(outputfolder,"RF_classif_repport_L2.txt")
# Define dataset to take into account
y_true=groundtruth_L2
y_pred=prediction_L2
class_label=classes_L2

# Compute precision accuracy
accuracy=accuracy_score(y_true, y_pred, normalize=True)
# Compute Cohen's Kappa
cohen_kappa=cohen_kappa_score(y_true, y_pred)
# Compute f1-score
f_1=f1_score(y_true, y_pred, average='weighted')
# Compute 'classification report'
classif_report=classification_report(y_true, y_pred, target_names=class_label)

# Save as .txt file
f=open(output, 'w')
f.write("Performance evalutation: \n")
f.write("Overall Accuracy: "+str(accuracy)+"\n")
f.write("Cohen's Kappa: "+str(cohen_kappa)+"\n")
f.write("F1-score: "+str(f_1)+"\n")
f.write("\n\n")
f.write("Classification report: \n "+classif_report)
f.close()

# Show file content
f=open(output,'r')
file_contents=f.read()
print (file_contents)
f.close()<jupyter_output>Performance evalutation: 
Overall Accuracy: 0.345
Cohen's Kappa: 0.272222222222
F1-score: 0.333011612458


Classification report: 
                             precision    recall  f1-score   support

                 Buildings       0.29      0.82      0.43        40
            Swimming pools       0.67      0.40      0.50        40
          Asphalt surfaces       0.58      0.38      0.45        40
       Brown/red bare soil       0.49      0.55      0.52        40
      White/grey bare soil       0.32      0.17      0.23        40
                     Trees       0.44      0.17      0.25        40
Mixed bare soil/vegetation       0.25      0.28      0.26        40
            Dry vegetation       0.19      0.10      0.13        40
          Other vegetation       0.20      0.33      0.25        40
             Inland waters       0.40      0.25      0.31        40

               avg / total       0.38      0.34      0.33       400<jupyter_text>Level 2b<jupyter_code>## Set the path to the output
output=os.path.join(outputfolder,"RF_classif_repport_L2b.txt")
# Define dataset to take into account
y_true=groundtruth_L2b
y_pred=prediction_L2b
class_label=classes_L2b

# Compute precision accuracy
accuracy=accuracy_score(y_true, y_pred, normalize=True)
# Compute Cohen's Kappa
cohen_kappa=cohen_kappa_score(y_true, y_pred)
# Compute f1-score
f_1=f1_score(y_true, y_pred, average='weighted')
# Compute 'classification report'
classif_report=classification_report(y_true, y_pred, target_names=class_label)

# Save as .txt file
f=open(output, 'w')
f.write("Performance evalutation: \n")
f.write("Overall Accuracy: "+str(accuracy)+"\n")
f.write("Cohen's Kappa: "+str(cohen_kappa)+"\n")
f.write("F1-score: "+str(f_1)+"\n")
f.write("\n\n")
f.write("Classification report: \n "+classif_report)
f.close()

# Show file content
f=open(output,'r')
file_contents=f.read()
print (file_contents)
f.close()<jupyter_output>Performance evalutation: 
Overall Accuracy: 0.455
Cohen's Kappa: 0.334351145038
F1-score: 0.44801237118


Classification report: 
                   precision    recall  f1-score   support

       Buildings       0.29      0.82      0.43        40
  Swimming pools       0.67      0.40      0.50        40
Asphalt surfaces       0.58      0.38      0.45        40
      Bare soils       0.49      0.41      0.45        80
  Low vegetation       0.53      0.57      0.55       120
           Trees       0.44      0.17      0.25        40
   Inland waters       0.40      0.25      0.31        40

     avg / total       0.49      0.46      0.45       400<jupyter_text>Level 1<jupyter_code>## Set the path to the output
output=os.path.join(outputfolder,"RF_classif_repport_L1.txt")
# Define dataset to take into account
y_true=groundtruth_L1
y_pred=prediction_L1
class_label=classes_L1

# Compute precision accuracy
accuracy=accuracy_score(y_true, y_pred, normalize=True)
# Compute Cohen's Kappa
cohen_kappa=cohen_kappa_score(y_true, y_pred)
# Compute f1-score
f_1=f1_score(y_true, y_pred, average='weighted')
# Compute 'classification report'
classif_report=classification_report(y_true, y_pred, target_names=class_label)

# Save as .txt file
f=open(output, 'w')
f.write("Performance evalutation: \n")
f.write("Overall Accuracy: "+str(accuracy)+"\n")
f.write("Cohen's Kappa: "+str(cohen_kappa)+"\n")
f.write("F1-score: "+str(f_1)+"\n")
f.write("\n\n")
f.write("Classification report: \n "+classif_report)
f.close()

# Show file content
f=open(output,'r')
file_contents=f.read()
print (file_contents)
f.close()<jupyter_output>Performance evalutation: 
Overall Accuracy: 0.5475
Cohen's Kappa: 0.347041847042
F1-score: 0.536448831


Classification report: 
                            precision    recall  f1-score   support

      Artificial surfaces       0.56      0.76      0.64       120
Natural material surfaces       0.49      0.41      0.45        80
               Vegetation       0.59      0.53      0.56       160
                    Water       0.40      0.25      0.31        40

              avg / total       0.54      0.55      0.54       400<jupyter_text>Reclass segmentation rasters with classification results<jupyter_code>## Print current mapset
print "You are currently working in the <"+str(mapsetname)+"> mapset"<jupyter_output>You are currently working in the <classif_C> mapset<jupyter_text>Change the outputfolder<jupyter_code>## Folder in which save processing time output
outputfolder=os.path.join(resultfolder,"classified_rasters")

## Create the folder if does not exists
if not os.path.exists(outputfolder):
    os.makedirs(outputfolder)
    print "Folder '"+outputfolder+"' created"<jupyter_output><empty_output><jupyter_text>Reclass segmentation raster Here after, the segment raster will be reclassed (using r.reclass) with the classification results.<jupyter_code># Connect to an existing database
db=pg.connect(database=dbname, user=dbuser, password=<PASSWORD>, host=host)
# Open a cursor to perform database operations
cur=db.cursor()
print "Start creating classified rasters"

classif_suffix=('l2','l2b','l1')
for x in classif_suffix:
    ## Set name for raster with the prediction at level 2
    prediction_raster='rf_classif_'+x

    #### Save 'seg_id' and the corresponding prediction for in a list
    query="SELECT seg_id, rf_pred_"+x+" FROM "+classifC_schema+"."+classif
    df=pd.read_sql(query, db)

    # Save seg_id of objects in a list
    listsegid=list(df['seg_id'])
    # Save predictions at level_2 in a list
    listpredict=list(df['rf_pred_'+x])

    if len(listsegid) <> len(listpredict):
        sys.exit('WARNING: lenght of lists containing segid and prediction are not the same. Please check before continue')
    else:
        nrows=len(listsegid)
        print str(nrows)+" segment are going to be reclassified"
    
    ##### Reclassify raster
    ## Create a temporary 'reclass_rule.csv' file
    temprulecsv=os.path.join(tempfile.gettempdir(),"reclass_rules.csv") # Define the csv output file name
    f = open(temprulecsv, 'w')
    # Write rules in the csv file
    for i in range(0,nrows-1):
        f.write(str(listsegid[i]))
        f.write("=")
        f.write(str(listpredict[i]))
        f.write("\n")
    f.write("*")
    f.write("=")
    f.write("NULL")
    f.close()

    ## Reclass segments raster layer to keep only outliers segments, using the reclas_rule.csv file (create temporary raster)
    grass.run_command('g.region', overwrite=True, raster="segments"+"@"+user["classificationA_mapsetname"])
    print ("Working on <"+prediction_raster+">: Reclassify original segment layer")
    grass.run_command('r.reclass', overwrite=True, input="segments"+"@"+user["classificationA_mapsetname"], 
                      output=prediction_raster, rules=temprulecsv)
    os.remove(temprulecsv)

    #### Create 'real raster'
    #### Make a copy of the classified maps of faster display in GRASS GIS
    ## Saving current time for processing time management
    print ("Working on <"+prediction_raster+">: Make hard copy of the reclassified layer")
    ## Set computational region
    grass.run_command('g.region', overwrite=True, raster="segments"+"@"+user["classificationA_mapsetname"])
    ## Create the same raster with r.mapcalc
    formula=prediction_raster+"_temp="+prediction_raster
    grass.mapcalc(formula, overwrite=True)
    ## Rename the new raster with the name of the original one (will be overwrited)
    print ("Working on <"+prediction_raster+">: Renaming layer")
    renameformula=prediction_raster+"_temp,"+prediction_raster
    grass.run_command('g.rename', overwrite=True, raster=renameformula)<jupyter_output>Start creating classified rasters
15438906 segment are going to be reclassified
Working on <rf_classif_l2>: Reclassify original segment layer
Working on <rf_classif_l2>: Make hard copy of the reclassified layer
Working on <rf_classif_l2>: Renaming layer
15438906 segment are going to be reclassified
Working on <rf_classif_l2b>: Reclassify original segment layer
Working on <rf_classif_l2b>: Make hard copy of the reclassified layer
Working on <rf_classif_l2b>: Renaming layer
15438906 segment are going to be reclassified
Working on <rf_classif_l1>: Reclassify original segment layer
Working on <rf_classif_l1>: Make hard copy of the reclassified layer
Working on <rf_classif_l1>: Renaming layer<jupyter_text>Change color table and export in .tif<jupyter_code># Define color table. Replace with the RGB values of wanted colors of each class
color_rule="11  227:26:28"+"\n"
color_rule+="12  255:141:1"+"\n"
color_rule+="13  94:221:227"+"\n"
color_rule+="14  102:102:102"+"\n"
color_rule+="21  246:194:142"+"\n"
color_rule+="22  211:217:173"+"\n"
color_rule+="31  0:128:0"+"\n"
color_rule+="32  189:255:185"+"\n"
color_rule+="33  88:190:141"+"\n"
color_rule+="34  29:220:0"+"\n"
color_rule+="41  30:30:192"+"\n"
color_rule+="51  0:0:0"+"\n"

## Create a temporary 'color_table.txt' file
color_table=os.path.join(outputfolder,"color_table_l2.txt") # Define the csv output file name
f = open(color_table, 'w')
f.write(color_rule)
f.close()
# Define color table. Replace with the RGB values of wanted colors of each class
color_rule="11  227:26:28"+"\n"
color_rule+="12  255:141:1"+"\n"
color_rule+="13  94:221:227"+"\n"
color_rule+="14  102:102:102"+"\n"
color_rule+="20  211:217:173"+"\n"
color_rule+="30  29:220:0"+"\n"
color_rule+="31  0:128:0"+"\n"
color_rule+="41  30:30:192"+"\n"
color_rule+="51  0:0:0"+"\n"

## Create a temporary 'color_table.txt' file
color_table=os.path.join(outputfolder,"color_table_l2b.txt") # Define the csv output file name
f = open(color_table, 'w')
f.write(color_rule)
f.close()
# Define color table. Replace with the RGB values of wanted colors of each class
color_rule="1   227:26:28"+"\n"
color_rule+="2   211:217:173"+"\n"
color_rule+="3   29:220:0"+"\n"
color_rule+="4   30:30:192"+"\n"
color_rule+="5   0:0:0"+"\n"

## Create a temporary 'color_table.txt' file
color_table=os.path.join(outputfolder,"color_table_l1.txt") # Define the csv output file name
f = open(color_table, 'w')
f.write(color_rule)
f.close()
# Create a list with classified raster
classifiedraster_list=grass.list_strings("rast", pattern="rf_classif_*", flag='r', mapset=classifC_schema)
print classifiedraster_list
for prediction_raster in classifiedraster_list:
    ## Apply new color the existing GRASS colortable (for faster display in GRASS map display)
    suffix=prediction_raster.split("@")[0].split("_")[-1]
    color_table_file=os.path.join(outputfolder,"color_table_"+suffix+".txt")
    grass.run_command('r.colors', map=prediction_raster, rules=color_table_file)<jupyter_output><empty_output><jupyter_text>Export of classification raster<jupyter_code>## Saving current time for processing time management
print ("Export classified raster maps on " + time.ctime())
begintime_exportraster=time.time()

for prediction_raster in classifiedraster_list:
    outputname=os.path.join(outputfolder,prediction_raster+".tif")
    grass.run_command('g.region', overwrite=True, raster=prediction_raster)
    grass.run_command('r.out.gdal', overwrite=True, input=prediction_raster, output=outputname, format='GTiff')
    
## Compute processing time and print it
print_processing_time(begintime_exportraster, "Classified raster maps exported in ")<jupyter_output>Export classified raster maps on Thu Dec  7 09:58:02 2017<jupyter_text>End of classification part<jupyter_code>print("The script ends at "+ time.ctime())
print_processing_time(begintime_segmentation_full, "Entire process has been achieved in ")<jupyter_output>The script ends at Mon Dec  4 16:53:01 2017
{'fraction_non_alphanumeric': 0.12279071258031951, 'fraction_numerical': 0.023442883009918758, 'mean_word_length': 5.172898120672602, 'pattern_counts': {'":': 0, '<': 252, '<?xml version=': 0, '>': 211, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '50', 'dir': 'jupyter-structured-clean-dedup', 'id': '214198', 'license': 'MIT', 'n_tokens_mistral': 21895, 'n_tokens_neox': 20709, 'n_words': 5290, 'path': 'Classification/Classif_C_SAR.ipynb', 'repo_name': 'tgrippa/Semi_automated_OBIA_processing_with_local_USPO', 'size': '62395'}