From 1d32fff14eab8a5a7edd5be5203f560090dfc199 Mon Sep 17 00:00:00 2001 From: Commandre Benjamin <benjamin.commandre@irstea.fr> Date: Fri, 7 Jun 2019 16:57:59 +0200 Subject: [PATCH] Regroupement des images en fonction de la tuile sentinel2 correspondante --- app/Archive.py | 105 ++++++++++++++++++---------------------------- app/Processing.py | 2 +- config.ini | 5 +++ main.py | 33 ++++++++------- 4 files changed, 64 insertions(+), 81 deletions(-) diff --git a/app/Archive.py b/app/Archive.py index 731aecc..98b0806 100644 --- a/app/Archive.py +++ b/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 diff --git a/app/Processing.py b/app/Processing.py index de521f0..c9dd1db 100644 --- a/app/Processing.py +++ b/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() diff --git a/config.ini b/config.ini index 32b6f8d..87b0415 100644 --- a/config.ini +++ b/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 diff --git a/main.py b/main.py index 07b3c50..d69c70d 100644 --- a/main.py +++ b/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): """ -- GitLab