From 632fa291ced2c89bd5d23683bc94540a50c36a69 Mon Sep 17 00:00:00 2001
From: Commandre Benjamin <benjamin.commandre@irstea.fr>
Date: Wed, 12 Sep 2018 17:07:13 +0200
Subject: [PATCH] Optimization + otb

---
 PHYMOBAT.py          |  5 ---
 Processing.py        | 58 +++++----------------------
 RasterSat_by_date.py |  5 ---
 Rpg.py               |  5 ---
 Sample.py            |  7 ----
 Segmentation.py      |  1 -
 Vector.py            | 95 +++++++++++---------------------------------
 7 files changed, 34 insertions(+), 142 deletions(-)

diff --git a/PHYMOBAT.py b/PHYMOBAT.py
index 9720371..992d507 100644
--- a/PHYMOBAT.py
+++ b/PHYMOBAT.py
@@ -77,11 +77,6 @@ class PHYMOBAT(QtWidgets.QMainWindow, Processing):
 			from ui_PHYMOBATs_tab import Ui_PHYMOBAT, _translate
 			
 		self.initUI()
-
-		# Chargement et lancement automatique pour les tests à supprimer
-		# self.open_backup("/media/commandre/8ce2dd88-1425-4161-a3c1-47f566c0285a/data_test/2018.xml")
-		# self.open_backup("/media/commandre/8ce2dd88-1425-4161-a3c1-47f566c0285a/Gironde_Test/gironde_total.xml")
-		# self.ok_button()
 		
 	def initUI(self):
 		
diff --git a/Processing.py b/Processing.py
index 96a2669..47fa560 100644
--- a/Processing.py
+++ b/Processing.py
@@ -43,6 +43,7 @@ from Vector import Vector
 from Sample import Sample
 import Segmentation
 from Rpg import Rpg
+import time
 
 import Constantes
 
@@ -299,28 +300,6 @@ class Processing(object):
 		for idx, i in enumerate(stats_name):
 			self.out_ndvistats_folder_tab[i] = path_stat[idx]
 
-		# self.logger.debug(self.out_ndvistats_folder_tab)
-		# raise(BaseException("Tst"))
-
-		# stats_ndvi, stats_cloud = current_list.calc_serie_stats(spectral_trans, self.folder_processing)
-
-		# # Stats cloud raster		
-		# stats_cloud_raster = "{0}/Cloud_number_{1}.TIF".format(self.folder_processing, self.captor_project)
-
-		# check_L8.complete_raster(\
-		# 	*check_L8.create_empty_raster(stats_cloud_raster, stats_cloud,\
-		# 	 check_L8.raster_data(self.check_download.list_img[0][4])[1]),\
-		# 	 stats_cloud)
-
-		# # Stats ndvi rasters		
-		# for stats_index in range(len(stats_ndvi)):
-		# 	out_ndvistats_raster = "{0}/{1}_{2}.TIF".format(self.folder_processing, stats_name[stats_index], self.captor_project)
-		# 	self.out_ndvistats_folder_tab[stats_index] = out_ndvistats_raster
-		# 	check_L8.complete_raster(\
-		# 		*check_L8.create_empty_raster(out_ndvistats_raster, stats_ndvi[stats_index],\
-		# 		 check_L8.raster_data(self.check_download.list_img[0][4])[1]),\
-		# 		 stats_ndvi[stats_index])
-
 		check_L8.logger.close()
 		current_list.logger.close()
 
@@ -547,6 +526,10 @@ class Processing(object):
 		self.logger.debug("Rasters : {0}".format(self.raster_path))
 		self.logger.debug(self.folder_processing)
 
+		self.logger.debug("sample verif : {0}".format(self.sample_name))
+
+		self.liste_chemin = [self.out_ndvistats_folder_tab['Min'], self.out_ndvistats_folder_tab['sfs'],\
+			self.out_ndvistats_folder_tab['haralick'], self.out_ndvistats_folder_tab['Max'] ]
 
 		# Extract value mean from polygons
 		for sple in range(len(self.sample_name) * 2):
@@ -575,22 +558,12 @@ class Processing(object):
 
 			items = list(self.raster_path.items())
 
-			self.logger.debug("Raster path : {0}".format(self.raster_path.items()))
+			self.logger.debug("Raster path items: {0}".format(self.raster_path.items()))
 
-			for lbo in range(len(self.raster_path)):
-				self.logger.debug("lbo : {0}".format(lbo))
-				kwargs['rank'] = lbo
-				kwargs['nb_img'] = len(self.raster_path)
-
-				self.logger.debug("item : {0}".format(items[lbo][1]))
-				sample_rd[sple].zonal_stats(items[lbo][1]["PATH"], items[lbo][1]["BAND"], **kwargs)
+			sample_rd[sple].zonal_stats(self.liste_chemin)
 
 			# To convert the dictionnary in a list
-			# for key, value in sample_rd[sple].stats_dict.iteritems():
 			for key, value in sample_rd[sple].stats_dict.items():
-
-				print("{0} : {1}".format(key, value))
-
 				# X_rf.append([-10000 if (math.isnan(x) or math.isinf(x)) else x for x in value])
 				# To set the grassland class of the RPG and PIAO (same class)			
 				if sple == 2 :
@@ -686,20 +659,9 @@ class Processing(object):
 		X_out_rf = [] # Variable list to compute decision tree with random forest method
 		if exist_stats == 0 :
 			out_carto.stats_dict = mgr.defaultdict(list)
-			
-			for i in range(len(multi_process_var)) :
-				kwargs['rank'] = i
-				kwargs['nb_img'] = len(multi_process_var)
-				p.append(Process(target=out_carto.zonal_stats, args=(*multi_process_var[i], ), kwargs=kwargs))
-				p[i].start()
-				
-				if self.mp == Constantes.MULTIPROCESSING_DISABLE:
-					p[i].join()
-			
-			if self.mp == Constantes.MULTIPROCESSING_ENABLE:	   
-				for i in range(len(multi_process_var)):
-					p[i].join()
-					
+
+			out_carto.zonal_stats(self.liste_chemin)
+								
 			for key, value_seg in out_carto.stats_dict.items():
 				true_value = [-10000 if (x is None or math.isnan(x) or math.isinf(x)) else x for x in value_seg]
 				X_out_rf.append(true_value)
diff --git a/RasterSat_by_date.py b/RasterSat_by_date.py
index d87702c..452dfcf 100644
--- a/RasterSat_by_date.py
+++ b/RasterSat_by_date.py
@@ -301,9 +301,6 @@ class RasterSat_by_date(object):
 
 		otb_MathBand = otb.Registry.CreateApplication("BandMathX")
 		otb_MathBand.SetParameterStringList("il", [self.cloudiness_pourcentage[0], self.cloudiness_pourcentage[1]])
-		# otb_MathBand.SetParameterString("exp",\
-		# "(im1b1 > -1000 ? (im2b1 == 0   ? (im1b{0}-im1b{1})/(im1b{0}+im1b{1}) : -1000) : im1b1)"\
-		# .format(canal_PIR, canal_R))
 		otb_MathBand.SetParameterString("exp",\
 		"(im1b1 > -1 ? (im2b1 == 0  ? ndvi(im1b{1}, im1b{0}) : -1) : im1b1)"\
 		.format(canal_PIR, canal_R))
@@ -332,8 +329,6 @@ class RasterSat_by_date(object):
 					  int -- variable **e**, Index to know if the raster exists.
 					  If it doesn't exists e = 0 else e = 1 (by default).
 		"""
-
-		# raise(BaseException("TEST"))
 		
 		e = 1 # Raster out exists by default 
 		
diff --git a/Rpg.py b/Rpg.py
index 7c09d25..95b858a 100644
--- a/Rpg.py
+++ b/Rpg.py
@@ -77,11 +77,6 @@ class Rpg(Vector):
         self.vector_used = "{0}/MONO_{1}".\
         format( os.path.dirname(self.vector_used), self.vector_name)
 
-        # self.logger.debug("self.vector_used : {0}".format(self.vector_used))
-        # self.logger.debug("self.vector_name : {0}".format(self.vector_name))
-        
-        # raise("TOP")
-
         if not os.path.exists(self.vector_used):
             
             shp_ogr = self.data_source.GetLayer()
diff --git a/Sample.py b/Sample.py
index b0ac80a..266127d 100644
--- a/Sample.py
+++ b/Sample.py
@@ -82,13 +82,6 @@ class Sample(Vector):
 		if not os.path.exists(self.output_folder) :
 			os.mkdir(self.output_folder)
 
-
-		"""
-			TODO : Enlever random.seed(0) pour retourner en aléatoire
-		"""
-
-		random.seed(0)
-
 		# If the users want select polygons with a certain class name
 		if kw_field and kw_classes:
 			# The random sample with class name selected only
diff --git a/Segmentation.py b/Segmentation.py
index 2ac1c11..0214a02 100644
--- a/Segmentation.py
+++ b/Segmentation.py
@@ -85,7 +85,6 @@ class Segmentation(Vector):
 		self.output_file =  'Final_classification.shp'
 		
 		self.out_class_name = []
-		# self.out_threshold = []
 		self.out_threshold = dict()
 
 		self.class_tab_final = defaultdict(list)
diff --git a/Vector.py b/Vector.py
index 8725f7c..3be93c3 100644
--- a/Vector.py
+++ b/Vector.py
@@ -134,57 +134,40 @@ class Vector():
 		
 		self.logger.info('Shapefile closing : {0}'.format(self.data_source.GetLayer().GetName()))
 	
-	def zonal_stats(self, inraster, band, **kwargs):
+	def zonal_stats(self, liste_chemin):
 		"""
-			Function to compute the mean in every polygons for a raster because of package ``rasterstats`` 
-			 in */usr/local/lib/python2.7/dist-packages/rasterstats-0.3.2-py2.7.egg/rasterstats/*
+			Function to compute the mean in every polygons on a list images with otb
 			
-			:param (inraster,band): inraster -> Input image path, and band -> band number
-			:type (inraster,band): tuple
-
-			:kwargs: **rank** (int) - Zonal stats ranking launch
-			
-					**nb_img** (int) - Number images launched with zonal stats 
+			:param liste_chemin : List input image path
+			:type liste_chemin: list(string)
 		"""
 		
-		
-		ranking = kwargs['rank'] if kwargs.get('rank') else 0
-		nb_img = kwargs['nb_img'] if kwargs.get('nb_img') else 1
-
-		time1 = time.time()
-		self.logger.info('Compute {0} stats on {1} band {2}'.format(os.path.split(self.vector_used)[1], os.path.split(inraster)[1], band))
+		self.logger.info("Compute stats 'mean' on {0}".format(os.path.split(self.vector_used)[1]))
 		
 		zonal_stats = otb.Registry.CreateApplication("ZonalStatistics")
-		zonal_stats.SetParameterString("in", inraster)
-		zonal_stats.SetParameterString("vec", self.vector_used)
-		zonal_stats.SetParameterInt("band", band-1)
-		zonal_stats.SetParameterInt("ram", 128)
-		statistique = ["mean"]
-		zonal_stats.SetParameterStringList("stats", statistique)
+		zonal_stats.SetParameterStringList("il", liste_chemin)
+		zonal_stats.SetParameterStringList("vec", [self.vector_used])
+		zonal_stats.SetParameterStringList("stats", ["mean"])
 
 		zonal_stats.Execute()
 
-		stats = list(zonal_stats.GetParameterValue("out"))
+		stats_otb = list(zonal_stats.GetParameterValue("out"))
 
-		# stats = rasterstats.zonal_stats(self.vector_used, inraster,  stats =['mean'], nodata=float('nan'), band=band) 
+		liste_valeur = []
 
-		for i in range(len(stats)):
-			temp = defaultdict(lambda : [0]*nb_img)
-			for j in range(nb_img):
-				try :
-					temp[0][j] = self.stats_dict[i][j]
-				except IndexError:
-					pass
+		for stats in stats_otb :
+			for s in stats.split(";"):
+				if s :
+					liste_valeur.append(np.array([float(x) for x in s.split()]))
 
-			# temp[0][ranking] = list(stats)[i]['mean']
-			temp[0][ranking] = float(stats[i])
-			self.stats_dict[i] = temp[0]
-			# self.logger.debug("self.stats_dict[i] : {0}".format(self.stats_dict[i]))
-		
-		time2 = time.time()
+		resultats = np.dstack(liste_valeur)[0]
+
+		stats_dict = defaultdict(list)
+
+		for idx, r in enumerate(resultats) :
+			self.stats_dict[idx] = r.tolist()
 					
-		self.logger.info('End of stats on {0} band {1}'.format(os.path.split(inraster)[1], band))
-		self.logger.debug("Temps d'éxecution {:.3f} ms".format((time2-time1)*1000.0))
+		self.logger.info("End of stats 'mean' on {0}".format(os.path.split(self.vector_used)[1]))
 
 	def zonal_stats_pp(self, inraster):
 		"""
@@ -195,41 +178,11 @@ class Vector():
 
 			:returns: dict -- **p_stats** : dictionnary with pxl percent in every polygon. Mainly 'Maj_count' (Majority Value) and 'Maj_count_perc' (Majority Percent)
 		"""
-		
-		# p_stats = rasterstats.zonal_stats(self.vector_used, inraster, stats=['count'], copy_properties=True, categorical=True)
-
-		# for i in range(len(p_stats)):
-		# 	percent = 0.0
-		# 	for p in p_stats[i].keys():
-		# 		if type(p) == np.float32 and p_stats[i][p]/float(p_stats[i]['count'])*100 > percent:
-		# 			p_stats[i]['Maj_count'] = p
-		# 			p_stats[i]['Maj_count_perc'] = p_stats[i][p]/float(p_stats[i]['count'])*100
-		# 			percent = p_stats[i]['Maj_count_perc']
-			
-		# 	if not 'Maj_count' in p_stats[i].keys() or not 'Maj_count_perc' in p_stats[i].keys():
-		# 		p_stats[i]['Maj_count']=0
-		# 		p_stats[i]['Maj_count_perc']=0
-
-		# for i in range(len(p_stats)):
-		# 	percent = 0.0
-		# 	Maj_count = 0
-		# 	Maj_count_perc = 0
-		# 	for p in p_stats[i].keys():
-		# 		if not isinstance(p, str) and  p_stats[i][p]/float(p_stats[i]['count'])*100 > percent:
-		# 			Maj_count = p
-		# 			Maj_count_perc = p_stats[i][p]/float(p_stats[i]['count'])*100
-		# 			percent = Maj_count_perc
-			
-		# 	p_stats[i]['Maj_count']= Maj_count
-		# 	p_stats[i]['Maj_count_perc']= Maj_count_perc
 
 		zonal_stats = otb.Registry.CreateApplication("ZonalStatistics")
-		zonal_stats.SetParameterString("in", inraster)
-		zonal_stats.SetParameterString("vec", self.vector_used)
-		zonal_stats.SetParameterInt("band", 1)
-		zonal_stats.SetParameterInt("ram", 128)
-		statistique = ["count"]
-		zonal_stats.SetParameterStringList("stats", statistique)
+		zonal_stats.SetParameterStringList("il", [inraster])
+		zonal_stats.SetParameterStringList("vec", [self.vector_used])
+		zonal_stats.SetParameterStringList("stats", ["count"])
 
 		zonal_stats.Execute()
 
-- 
GitLab