Amélioration du code

app/Archive.py

@@ -15,6 +15,7 @@ from uuid import uuid4
 from collections import defaultdict, UserDict
 import app.Satellites as Satellites
 import app.Outils as Outils
+import app.Constantes as Constantes
 
 class Archive():
     """
@@ -39,14 +40,23 @@ class Archive():
         :type annee_fin: Entier
     """
 
-    def __init__(self, capteur, extensions, niveau, emprise, sortie, annee_debut, annee_fin):
+    def __init__(self, capteur, bandes, niveau, emprise, zone_etude, sortie, annee_debut, annee_fin, seuil):
         """
             Créé une instance de la classe 'Archive'.
         """
 
         self._capteur = capteur
         self.niveau = niveau
-        self.extent_img = extensions
+       
+        if bandes == "RGB" :
+            self.extent_img = Satellites.SATELLITE[self._capteur][niveau][:-1]
+        else :
+            self.extent_img = Satellites.SATELLITE[self._capteur][niveau]
+
+        if self.niveau == "LEVEL2A" :
+            self.nuage = Satellites.SATELLITE[self._capteur]["NUAGE"]
+        else : 
+            self.nuage = None
 
         self.liste_annees = []
 
@@ -56,9 +66,10 @@ class Archive():
         self.liste_annees.append(annee_fin)
 
         self.emprise = emprise
+        self.zone_etude = zone_etude
         self.dossier_sortie = sortie
-        self.emprise = emprise
-    
+        self.seuil_nuage = seuil
+
         self.liste_archive = []
 
         self.logger = Outils.Log("log", "Archive")
@@ -392,6 +403,43 @@ class Archive():
 
         self.logger.info("All images have been downloaded !")
 
+    def calcul_ennuagement(self, tuiles_nuage):
+
+        option_warp     = gdal.WarpOptions(format='Mem', cropToCutline=True, outputType=gdal.GDT_Float32,\
+        								   cutlineDSName=self.emprise, dstNodata="NaN")
+        mosaic_nuage    = gdal.Warp("", tuiles_nuage, options=option_warp)
+
+        rows = mosaic_nuage.RasterYSize # Rows number
+        cols = mosaic_nuage.RasterXSize # Columns number
+
+        data = mosaic_nuage.GetRasterBand(1).ReadAsArray(0, 0, cols, rows)
+
+        nombre_pixel_zone = 0
+        nombre_pixel_null = 0
+        nombre_pixel_non_nuage = 0
+
+        for donnees in data :
+            d = donnees.tolist()
+            nombre_pixel_zone += len(d)
+            nombre_pixel_null += d.count(-10000) + np.isnan(d).sum()
+            nombre_pixel_non_nuage += d.count(0)
+
+        nombre_pixel_non_null = nombre_pixel_zone - nombre_pixel_null
+
+        # self.logger.debug("nombre_pixel_non_null : {}".format(nombre_pixel_non_null))
+        # self.logger.debug("nombre_pixel_null : {}".format(nombre_pixel_null))
+        # self.logger.debug("nombre_pixel_non_nuage : {}".format(nombre_pixel_non_nuage))
+
+        del data
+
+        ennuagement = (float(nombre_pixel_non_nuage)/float(nombre_pixel_non_null)) if nombre_pixel_non_null != 0 else 0.0
+        self.logger.debug("Ennuagement : {}%".format(100 - round(ennuagement*100, 2)))
+   
+        # Computer cloud's percentage with dist (sum of cloud) by sum of the image's extent
+        # return mosaic_nuage, 1.0 - (float(dist)/(np.sum(mask_spec)) if np.sum(mask_spec) != 0 else 0)
+        return 1.0 - ennuagement
+
+
     def calcul_aire(self, tuiles_image):
         """
             Méthode pour fusionner un ensemble d'images et calculer l'aire de l'image obtenue.
@@ -409,10 +457,10 @@ class Archive():
         cols = mosaic_image.RasterXSize # Columns number
 
         data        = mosaic_image.GetRasterBand(1).ReadAsArray(0, 0, cols, rows).astype(np.float16)
-        mask_spec   = np.isin(data, [-10000, math.isnan], invert=True)
+        mask_spec   = ~(np.isnan(data) | np.isin(data, [-10000]))
 
         area = float((np.sum(mask_spec) * mosaic_image.GetGeoTransform()[1] * abs(mosaic_image.GetGeoTransform()[-1]) )/10000)
-        self.logger.info("Aire = {0} ha".format(area))
+        # self.logger.info("Aire = {0} ha".format(area))
 
         return mosaic_image, area
 
@@ -458,10 +506,11 @@ class Archive():
         self.logger.info("Date : {0} -> {1} image(s)".format(date, len(liste_content)))
 
         tuiles_image = []
+        tuiles_nuage = []
 
         # Options Gdal pour la fusion des bandes
-        options_vrt         = gdal.BuildVRTOptions(separate=True)
-        options_translate   = gdal.TranslateOptions(format="Mem")
+        options_vrt     	= gdal.ParseCommandLine('-resolution highest -srcnodata -10000 -vrtnodata NaN -separate')
+        options_translate   = gdal.ParseCommandLine('-of Mem -ot Float32 -a_nodata NaN')
 
         self.logger.info("Extraction des images")
 
@@ -493,12 +542,32 @@ class Archive():
             # on découpe l'image selon l'emprise
             # et on conserve l'image obtenue
             vrt = gdal.BuildVRT("", liste_bandes, options=options_vrt)
-            tuiles_image.append(Outils.clip(gdal.Translate("", vrt, options=options_translate), self.emprise))
 
-            # On libère la mémoire
+            image = Outils.clip(vrt, self.zone_etude)
+            _, aire = self.calcul_aire(image)
+
+            if aire > 0.0 :
+                tuiles_image.append(vrt)
+
+                # On libère la mémoire
+                for mmap_name in liste_mem :
+                    gdal.Unlink(mmap_name)
+
+                liste_mem = []
+
+                if self.nuage is not None :
+                    image_nuage = [f for f in zip_img if self.nuage[0] in f][0]
+                    mmap_name = "/vsimem/"+uuid4().hex
+
+                    gdal.FileFromMemBuffer(mmap_name, tzip.read(image_nuage))
+                    liste_mem.append(mmap_name)
+                    tuiles_nuage.append(Outils.clip(gdal.Open(mmap_name), self.zone_etude))
+
             for mmap_name in liste_mem :
                 gdal.Unlink(mmap_name)
 
+            liste_mem = []
+
             liste_bandes = None
 
             del tzip
@@ -506,21 +575,22 @@ class Archive():
 
         del liste_content
 
-        # On effectue une mosaïque des images qu'on découpe selon l'emprise,
-        # puis on calcule l'aire de l'image ainsi obtenue
-        self.logger.info("Calcul de l'aire")
-        image, aire = self.calcul_aire(tuiles_image)
-        del tuiles_image
+        if len(tuiles_image) > 0 :          
+
+            if len(tuiles_nuage) == 0 or self.calcul_ennuagement(tuiles_nuage) <= self.seuil_nuage:
+
+                # On effectue une mosaïque des images qu'on découpe selon l'emprise.
+                image = Outils.clip(tuiles_image, self.emprise)
 
-        # Si l'aire est positive on enregistre l'image
-        if aire > 0.0 :
+                self.logger.info("Sauvegarde des images")
 
-            self.logger.info("Sauvegarde des images")
+                dossier = "{0}/{1}/Images".format(self.dossier_sortie, date[:4])
 
-            dossier = "{0}/{1}/Images".format(self.dossier_sortie, date[:4])
+                self.logger.debug("Dossier image : {0}".format(dossier))
 
-            self.logger.debug("Dossier image : {0}".format(dossier))
+                self.ecriture_geotiff(image, "{0}/{1}.tif".format(dossier,date))
 
-            self.ecriture_geotiff(image, "{0}/{1}.tif".format(dossier,date))
+                del tuiles_nuage
 
-        del image
\ No newline at end of file
+            del tuiles_image
+            del image
\ No newline at end of file

app/Constantes.py

@@ -4,13 +4,4 @@
 import logging
 
 NIVEAU_DEFAUT 	= logging.DEBUG
-
-CLOUD_THRESHOLD = 0.0
-
 EPSG_PHYMOBAT 	= 2154
-
-# Extensions des bandes bleue, verte et rouge
-BANDES_RGB = ['_FRC_B2.tif', '_FRC_B3.tif', '_FRC_B4.tif']
-
-# Extensions des bandes bleue, verte et rouge + proche infra-rouge
-BANDES_RGBI = ['_FRC_B2.tif', '_FRC_B3.tif', '_FRC_B4.tif', '_FRC_B8.tif']

app/Outils.py

@@ -97,7 +97,7 @@ def limitation_temporelle(secondes):
 # except TimeoutException:
 #     pass
 
-def clip(image, cut, form="Mem", dst=""):
+def clip(image, cut, form="Mem", dst="", type=gdal.GDT_Float32, nodata="NaN"):
 	"""
 		Découpe une image selon un shapefile donné. Si le type de l'image en sortie
 		n'est pas spécifié, l'image sera retournée en pointeur mémoire et ne sera pas sauvegardée.
@@ -116,6 +116,6 @@ def clip(image, cut, form="Mem", dst=""):
 	"""
         
 	option_clip = gdal.WarpOptions(cropToCutline=True,\
-	cutlineDSName=cut, outputType=gdal.GDT_Float32 , format=form,  dstNodata='NaN')
+	cutlineDSName=cut, outputType=type , format=form, dstNodata=nodata)
 
 	return gdal.Warp(dst, image, options=option_clip)
\ No newline at end of file

app/Processing.py

@@ -27,7 +27,8 @@ class Processing(object):
 		if not self.annee_fin :
 			self.annee_fin = datetime.datetime.now().year
 
-		self.check_download = Archive.Archive(self.capteur, self.extensions, self.niveau, self.emprise, self.resultats, self.annee_debut, int(self.annee_fin))
+		self.check_download = Archive.Archive(self.capteur, self.bandes, self.niveau, self.emprise, self.zone_etude,\
+											  self.resultats, self.annee_debut, int(self.annee_fin), self.seuil_nuage)
 		self.check_download.listing()
 		self.check_download.download_auto(self.id, self.mdp, self.proxy)
 
@@ -99,5 +100,5 @@ class Processing(object):
 		"""
 			Calul le ndvi, fusionnne en une seule image puis lance le module OTBPhenology
 		"""
-		self.calcul_ndvi()
-		self.otbPhenologie()
\ No newline at end of file
+		# self.calcul_ndvi()
+		# self.otbPhenologie()
\ No newline at end of file

app/Satellites.py

@@ -14,6 +14,9 @@ SATELLITE["SENTINEL2"]["resto"] = "resto2"
 SATELLITE["SENTINEL2"]["token_type"] = "text"
 SATELLITE["SENTINEL2"]["R"] = 2
 SATELLITE["SENTINEL2"]["PIR"] = 3
+SATELLITE["SENTINEL2"]["LEVEL2A"] 	= ['_SRE_B2.tif', '_SRE_B3.tif', '_SRE_B4.tif', '_SRE_B8.tif']
+SATELLITE["SENTINEL2"]["NUAGE"] 	= ['_CLM_R1.tif']
+SATELLITE["SENTINEL2"]["LEVEL3A"] 	= ['_FRC_B2.tif', '_FRC_B3.tif', '_FRC_B4.tif', '_FRC_B8.tif']
 
 ########################## LANDSAT #################################
 

main.py

@@ -34,11 +34,7 @@ class Telechargement(Processing):
 		# Capteur utilisé
 		self.capteur = "{}".format(configfile["satellite"]["capteur"])
 		self.niveau = "{}".format(configfile["satellite"]["processingLevel"])
-
-		if configfile["satellite"]["bandes"] == "RGB" :
-			self.extensions = Constantes.BANDES_RGB
-		else :
-			self.extensions = Constantes.BANDES_RGBI
+		self.bandes = "{}".format(configfile["satellite"]["bandes"])
 		
 		# Date de début et de fin de la recherche
 		try:
@@ -47,9 +43,14 @@ class Telechargement(Processing):
 			raise "L'année de départ est requise."
 		
 		self.annee_fin = "{}".format(configfile["donnees"]["annee_fin"])
+		self.seuil_nuage = float("{}".format(configfile["donnees"]["seuil_nuage"]))/100.0
 		
-		# Emprise de l'étude
-		self.emprise = "{}".format(configfile["donnees"]["chemin_emprise"])
+		# Emprise et zone de l'étude
+		self.emprise 	= "{}".format(configfile["donnees"]["chemin_emprise"])
+		self.zone_etude = "{}".format(configfile["donnees"]["chemin_zone_etude"])
+
+		if not self.zone_etude :
+			self.zone_etude = self.emprise
 		
 		# Identifiant, mot de passe et proxy pour le téléchargement des images Théia
 		self.id 	= "{}".format(configfile["theia"]["identifiant"])
@@ -60,9 +61,8 @@ class Telechargement(Processing):
 	def run(self):
 		"""
 			Fonction pour lancer le programme
-		"""
-		
-		# Début du processus
+		"""		
+		Début du processus
 		debut = time.time()
 		
 		# Recherche de nouvelles images non traitées et téléchargement de celles-ci le cas échéant
@@ -77,7 +77,16 @@ class Telechargement(Processing):
 		self.logger.info('Programme terminé en {} secondes'.format(fin - debut))
 		nb_jours = int(time.strftime('%d', time.gmtime(fin - debut))) - 1
 		self.logger.info("Cela représente {} jour(s) {}".format(nb_jours, time.strftime('%Hh %Mmin%S', time.gmtime(fin - debut))))
-		
+
 if __name__ == "__main__":
+
 	app = Telechargement()
-	sys.exit(app.run())
\ No newline at end of file
+
+	SUCCESS = False
+
+	while not SUCCESS :
+		try:
+			sys.exit(app.run())
+			SUCCESS = True
+		except Exception as e:
+			pass	
\ No newline at end of file