Regroupement des images en fonction de la tuile sentinel2 correspondante

app/Archive.py

@@ -8,7 +8,7 @@ import zipfile
 import requests
 import configparser
 
-from osgeo import ogr, gdal
+from osgeo import gdal, ogr, osr
 import numpy as np
 from uuid import uuid4
 
@@ -17,6 +17,9 @@ import app.Satellites as Satellites
 import app.Outils as Outils
 import app.Constantes as Constantes
 
+id_tuile   = re.compile("T[0-9]+[A-Z]+")
+date_tuile = re.compile("[SENTINEL.+?_]([0-9]*?)-")
+
 class Archive():
     """
         Classe pour lister, télécharger les archives via site Theia selon un shapefile représentant la zone d'étude.
@@ -24,12 +27,18 @@ class Archive():
         :param capteur: Nom du satellite utilisé (ex: Landsat, Sentinel2, ...).
         :type capteur: Chaîne de caractères
 
+        :param bandes: Bandes à conservées.
+        :type bandes: Chaîne de caractères
+
         :param niveau: Niveau de traitement voulu.
         :type niveau: Chaîne de caractères
 
-        :param emprise: Chemin vers le shapefile représentant la zone d'étude.
+        :param emprise: Chemin vers le shapefile représentant l'emprise.
         :type emprise: Chaîne de caractères
 
+        :param zone_etude: Chemin vers le shapefile représentant la zone d'étude.
+        :type zone_etude: Chaîne de caractères
+
         :param sortie: Chemin du dossier de sortie.
         :type sortie: Chaîne de caractères
 
@@ -38,6 +47,9 @@ class Archive():
 
         :param annee_fin: Année jusqu'à laquelle on cherche les images.
         :type annee_fin: Entier
+
+        :param seuil: Seuil maximal d'ennuagement.
+        :type seuil: Entier
     """
 
     def __init__(self, capteur, bandes, niveau, emprise, zone_etude, sortie, annee_debut, annee_fin, seuil):
@@ -61,7 +73,7 @@ class Archive():
         self.liste_annees = []
 
         for i in range(annee_debut, annee_fin) :
-             self.liste_annees.append(i)
+            self.liste_annees.append(i)
 
         self.liste_annees.append(annee_fin)
 
@@ -80,7 +92,7 @@ class Archive():
 
         self.url = ''
 
-    def utm_to_latlng(self, zone, easting, northing, northernHemisphere=True):
+    def utm_to_latlng(self, zone, easting, northing):
         """
             Méthode pour convertir la projection UTM en coordonnées géographique
 
@@ -97,57 +109,23 @@ class Archive():
             :type northernHemisphere: Booléen
 
             :returns: Tuple d'entiers représentant la longitude et la latitude.
-
-            Source : http://www.ibm.com/developerworks/java/library/j-coordconvert/index.html
         """
 
-        if not northernHemisphere:
-            northing = 10000000 - northing
-
-        a = 6378137
-        e = 0.081819191
-        e1sq = 0.006739497
-        k0 = 0.9996
-
-        arc = northing / k0
-        mu = arc / (a * (1 - math.pow(e, 2) / 4.0 - 3 * math.pow(e, 4) / 64.0 - 5 * math.pow(e, 6) / 256.0))
-
-        ei = (1 - math.pow((1 - e * e), (1 / 2.0))) / (1 + math.pow((1 - e * e), (1 / 2.0)))
-
-        ca = 3 * ei / 2 - 27 * math.pow(ei, 3) / 32.0
-
-        cb = 21 * math.pow(ei, 2) / 16 - 55 * math.pow(ei, 4) / 32
-        cc = 151 * math.pow(ei, 3) / 96
-        cd = 1097 * math.pow(ei, 4) / 512
-        phi1 = mu + ca * math.sin(2 * mu) + cb * math.sin(4 * mu) + cc * math.sin(6 * mu) + cd * math.sin(8 * mu)
-
-        n0 = a / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (1 / 2.0))
+        systeme_utm = osr.SpatialReference()
 
-        r0 = a * (1 - e * e) / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (3 / 2.0))
-        fact1 = n0 * math.tan(phi1) / r0
+        # Initialise le système de coordonnées géographique vers Lat/Long
+        systeme_utm.SetWellKnownGeogCS("WGS84")
+        systeme_utm.SetUTM(zone, northing > 0)
 
-        _a1 = 500000 - easting
-        dd0 = _a1 / (n0 * k0)
-        fact2 = dd0 * dd0 / 2
+        # Duplique UNIQUEMENT le système de coordonnées géographique
+        systeme_lat_long = systeme_utm.CloneGeogCS()
 
-        t0 = math.pow(math.tan(phi1), 2)
-        Q0 = e1sq * math.pow(math.cos(phi1), 2)
-        fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * math.pow(dd0, 4) / 24
+        # Creation de méthode de conversion
+        # (<de>, <vers>)
+        convertion_utm = osr.CoordinateTransformation(systeme_utm, systeme_lat_long)
 
-        fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * math.pow(dd0, 6) / 720
-
-        lof1 = _a1 / (n0 * k0)
-        lof2 = (1 + 2 * t0 + Q0) * math.pow(dd0, 3) / 6.0
-        lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * math.pow(Q0, 2) + 8 * e1sq + 24 * math.pow(t0, 2)) * math.pow(dd0, 5) / 120
-        _a2 = (lof1 - lof2 + lof3) / math.cos(phi1)
-        _a3 = _a2 * 180 / math.pi
-
-        latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / math.pi
-
-        if not northernHemisphere:
-            latitude = -latitude
-
-        longitude = ((zone > 0) and (6 * zone - 183.0) or 3.0) - _a3
+        # Renvoie lon, lat, altitude
+        longitude, latitude, altitude = convertion_utm.TransformPoint(easting, northing, 0)
 
         return (longitude, latitude)
 
@@ -187,10 +165,10 @@ class Archive():
         data_source.Destroy()
 
         # Coordinates min in decimal degrees
-        LL_min = self.utm_to_latlng(31, extent_[0],  extent_[2], northernHemisphere=True)
+        LL_min = self.utm_to_latlng(31, extent_[0],  extent_[2])
 
         # Coordinates max in decimal degrees
-        LL_max = self.utm_to_latlng(31, extent_[1],  extent_[3], northernHemisphere=True)
+        LL_max = self.utm_to_latlng(31, extent_[1],  extent_[3])
 
         area_coord_corner = str(LL_min[0]) + ',' + str(LL_min[1]) + ',' + str(LL_max[0]) + ',' + str(LL_max[1])
 
@@ -261,7 +239,7 @@ class Archive():
 
         self.logger.info("{0} image(s) correspondent aux critères.".format(len(self.liste_archive)))
 
-    def download_auto(self, identifiant, mdp, proxy="", extraction=True):
+    def download_auto(self, identifiant, mdp, proxy="", extraction=True, groupe="Date"):
         """
             Méthode pour télécharger les archives sur le site Theia.
 
@@ -308,8 +286,6 @@ class Archive():
 
         head = {"Authorization": "Bearer {0}".format(token)}
 
-        regex = 'SENTINEL(.+?)_(.+?)-'
-
         # Dictionnaire des images à traiter regrouper par date
         dico_date = dict()
 
@@ -319,7 +295,7 @@ class Archive():
         for l in self.liste_archive :
 
             # Récupération de la date quand l'image a été prise
-            date = re.search(regex, l[1]).group(2)
+            date = date_tuile.search(l[1]).group(1)
 
             # Récupération du nom de l'image
             nom_image = l[1].split("/")[-1][:-4]
@@ -386,8 +362,16 @@ class Archive():
                     liste_content.append(reponse.content)
                     del reponse
                 else :
-                    dossier = "{0}/{1}/Images".format(self.dossier_sortie, cle[:4])
-                    with open("{0}/{1}".format(dossier,img[1].split("/")[-1]), "wb") as fichier:
+                    if groupe == "Tuile" :
+                        dossier_archive = "{0}/{1}/Archive".format(self.dossier_sortie, id_tuile.search(img[1]).group(0))
+                    else:
+                        dossier_archive = "{0}/{1}/Archive".format(self.dossier_sortie, cle[:4])
+
+
+                    if not os.path.exists(dossier_archive):
+                        os.makedirs(dossier_archive)
+
+                    with open("{0}/{1}".format(dossier_archive, img[1].split("/")[-1]), "wb") as fichier:
                         fichier.write(reponse.content)
 
             if extraction :
@@ -433,20 +417,14 @@ class Archive():
 
         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.
@@ -467,7 +445,6 @@ class Archive():
         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))
 
         return mosaic_image, area
 

app/Processing.py

@@ -30,7 +30,7 @@ class Processing(object):
 		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, extraction=self.extraction)
+		self.check_download.download_auto(self.id, self.mdp, self.proxy, extraction=self.extraction, groupe=self.groupe)
 
 		self.liste_dossier = dict()
 

config.ini

@@ -20,6 +20,11 @@ seuil_nuage = 5.0
 # chemin/dossier/resultats
 chemin =
 
+# Manière de regrouper les images/archives
+# Possibilités :- Date (Regroupe les images selon leurs date d'acquisition)
+#		- Tuile (Regroupe les images en fonction de la tuile sentinel2 auquelle elles correspondent )
+groupe = Date
+
 # Vrai si extraction des images des archives, faux si téléchagement de celles-ci
 extraction = True
 

main.py

@@ -31,43 +31,44 @@ class Telechargement(Processing):
 		configfile = configparser.ConfigParser()
 		configfile.read("config.ini")
 
-		# Dossier contenant les résultats
-		self.resultats = "{}".format(configfile["sortie"]["chemin"])
-
 		# Capteur utilisé
-		self.capteur = "{}".format(configfile["satellite"]["capteur"])
-		self.niveau = "{}".format(configfile["satellite"]["processingLevel"])
-		self.bandes = "{}".format(configfile["satellite"]["bandes"])
+		self.capteur = configfile["satellite"]["capteur"]
+		self.niveau  = configfile["satellite"]["processingLevel"]
+		self.bandes  = configfile["satellite"]["bandes"]
+
+		# Dossier contenant les résultats
+		self.resultats = configfile["sortie"]["chemin"]
 
 		try:
-			if str2bool("{}".format(configfile["sortie"]["extraction"])):
+			if str2bool(configfile["sortie"]["extraction"]):
 				self.extraction = False
 			else :
 				self.extraction = True
 		except :
 			self.extraction = True
 
+		self.groupe = configfile["sortie"]["groupe"]
+
 		# Date de début et de fin de la recherche
 		try:
-			self.annee_debut = int("{}".format(configfile["donnees"]["annee_debut"]))
+			self.annee_debut = int(configfile["donnees"]["annee_debut"])
 		except Exception as e:
 			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 if "{}".format(configfile["donnees"]["seuil_nuage"]) else 0.0
+		self.annee_fin = configfile["donnees"]["annee_fin"]
+		self.seuil_nuage = float(configfile["donnees"]["seuil_nuage"])/100.0 if configfile["donnees"]["seuil_nuage"] else 0.0
 
 		# Emprise et zone de l'étude
-		self.emprise 	= "{}".format(configfile["donnees"]["chemin_emprise"])
-		self.zone_etude = "{}".format(configfile["donnees"]["chemin_zone_etude"])
+		self.emprise 	= configfile["donnees"]["chemin_emprise"]
+		self.zone_etude = 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"])
-		self.mdp 	= "{}".format(configfile["theia"]["mdp"])
-		self.proxy 	= "{}".format(configfile["theia"]["proxy"])
-
+		self.id 	= configfile["theia"]["identifiant"]
+		self.mdp 	= configfile["theia"]["mdp"]
+		self.proxy 	= configfile["theia"]["proxy"]
 
 	def run(self):
 		"""