Commit 03ea8e75 authored by Commandre Benjamin's avatar Commandre Benjamin
Browse files

Optimisation + correction du calcul des coordonnees

parent 1d32fff1
...@@ -6,20 +6,24 @@ import math, subprocess ...@@ -6,20 +6,24 @@ import math, subprocess
import urllib.request import urllib.request
import zipfile import zipfile
import requests import requests
import datetime
import configparser import configparser
from urllib.parse import urlencode
from osgeo import gdal, ogr, osr from osgeo import gdal, ogr, osr
import numpy as np import numpy as np
from uuid import uuid4 from uuid import uuid4
from collections import defaultdict, UserDict from collections import defaultdict, UserDict
import app.Satellites as Satellites
import app.Outils as Outils import app.Outils as Outils
import app.Constantes as Constantes import app.Constantes as Constantes
id_tuile = re.compile("T[0-9]+[A-Z]+") id_tuile = re.compile("T[0-9]+[A-Z]+")
date_tuile = re.compile("[SENTINEL.+?_]([0-9]*?)-") date_tuile = re.compile("[SENTINEL.+?_]([0-9]*?)-")
def str2bool(v):
return v.lower() in (["false"])
class Archive(): class Archive():
""" """
Classe pour lister, télécharger les archives via site Theia selon un shapefile représentant la zone d'étude. Classe pour lister, télécharger les archives via site Theia selon un shapefile représentant la zone d'étude.
...@@ -52,82 +56,70 @@ class Archive(): ...@@ -52,82 +56,70 @@ class Archive():
:type seuil: Entier :type seuil: Entier
""" """
def __init__(self, capteur, bandes, niveau, emprise, zone_etude, sortie, annee_debut, annee_fin, seuil): # def __init__(self, capteur, bandes, niveau, emprise, zone_etude, sortie, annee_debut, annee_fin, seuil):
def __init__(self, config):
""" """
Créé une instance de la classe 'Archive'. Créé une instance de la classe 'Archive'.
""" """
self._capteur = capteur
self.niveau = niveau
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 = []
for i in range(annee_debut, annee_fin) :
self.liste_annees.append(i)
self.liste_annees.append(annee_fin)
self.emprise = emprise
self.zone_etude = zone_etude
self.dossier_sortie = sortie
self.seuil_nuage = seuil
self.liste_archive = []
self.logger = Outils.Log("log", "Archive") self.logger = Outils.Log("log", "Archive")
self.serveur = Satellites.SATELLITE[self._capteur]["serveur"] Outils.get_variable(self, config)
self.resto = Satellites.SATELLITE[self._capteur]["resto"]
self.token_type = Satellites.SATELLITE[self._capteur]["token_type"]
self.liste_archive = []
self.url = '' self.url = ''
def utm_to_latlng(self, zone, easting, northing): def vector_projection_with_epsg(self, chemin_vecteur, output_path, epsg=32631):
"""
Méthode pour convertir la projection UTM en coordonnées géographique driver = ogr.GetDriverByName('ESRI Shapefile')
:param zone: Zone en projection UTM. projection_reference = osr.SpatialReference()
:type zone: Entier projection_reference.ImportFromEPSG(epsg)
:param easting: Coordonnées UTM en x. shapefile = driver.Open(chemin_vecteur)
:type easting: Flottant layer = shapefile.GetLayer()
projection_vector = layer.GetSpatialRef()
:param northing: Coordonnées UTM en y. transformation = osr.CoordinateTransformation(projection_vector, projection_reference)
:type northing: Flottant
:param northernHemisphere: Vrai si la zone se situe dans l'hémisphère Nord, faux sinon. vecteur_reproj = output_path
:type northernHemisphere: Booléen
:returns: Tuple d'entiers représentant la longitude et la latitude. if os.path.exists(vecteur_reproj):
""" driver.DeleteDataSource(vecteur_reproj)
systeme_utm = osr.SpatialReference() outDataSet = driver.CreateDataSource(vecteur_reproj)
outLayer = outDataSet.CreateLayer(layer.GetName(), projection_reference, geom_type=layer.GetLayerDefn().GetGeomType())
# Initialise le système de coordonnées géographique vers Lat/Long # add fields
systeme_utm.SetWellKnownGeogCS("WGS84") inLayerDefn = layer.GetLayerDefn()
systeme_utm.SetUTM(zone, northing > 0) for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
outLayer.CreateField(fieldDefn)
# Duplique UNIQUEMENT le système de coordonnées géographique outLayerDefn = outLayer.GetLayerDefn()
systeme_lat_long = systeme_utm.CloneGeogCS()
# Creation de méthode de conversion inFeature = layer.GetNextFeature()
# (<de>, <vers>)
convertion_utm = osr.CoordinateTransformation(systeme_utm, systeme_lat_long)
# Renvoie lon, lat, altitude while inFeature:
longitude, latitude, altitude = convertion_utm.TransformPoint(easting, northing, 0) # get the input geometry
geom = inFeature.GetGeometryRef()
# reproject the geometry
geom.Transform(transformation)
# create a new feature
outFeature = ogr.Feature(outLayerDefn)
# set the geometry and attribute
outFeature.SetGeometry(geom)
for i in range(0, outLayerDefn.GetFieldCount()):
outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
# add the feature to the shapefile
outLayer.CreateFeature(outFeature)
# dereference the features and get the next input feature
outFeature = None
inFeature = layer.GetNextFeature()
return (longitude, latitude) # Save and close the shapefiles
inDataSet = None
outDataSet = None
def coord_box_dd(self): def coord_box_dd(self):
""" """
...@@ -144,8 +136,8 @@ class Archive(): ...@@ -144,8 +136,8 @@ class Archive():
if os.path.exists(utm_outfile): if os.path.exists(utm_outfile):
os.remove(utm_outfile) os.remove(utm_outfile)
process_tocall = ['ogr2ogr', '-f', 'Esri Shapefile', '-t_srs', 'EPSG:32631', utm_outfile, self.emprise] # wgs84 : 4326
subprocess.call(process_tocall) self.vector_projection_with_epsg(self.emprise, utm_outfile, epsg=4326)
# To get shapefile extent # To get shapefile extent
# Avec import ogr # Avec import ogr
...@@ -164,15 +156,17 @@ class Archive(): ...@@ -164,15 +156,17 @@ class Archive():
## Close source data ## Close source data
data_source.Destroy() data_source.Destroy()
# Coordinates min in decimal degrees area_coord_corner = str(extent_[2]) + ',' + str(extent_[0]) + ',' + str(extent_[3]) + ',' + str(extent_[1])
LL_min = self.utm_to_latlng(31, extent_[0], extent_[2])
return "{},{},{},{}".format(extent_[2], extent_[0], extent_[3], extent_[1])
def listing_by_tile(self):
# Coordinates max in decimal degrees for tuile in self.liste_tuiles :
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]) self.requete['location'] = tuile
return area_coord_corner print(urlencode(self.requete))
def listing(self): def listing(self):
""" """
...@@ -181,15 +175,33 @@ class Archive(): ...@@ -181,15 +175,33 @@ class Archive():
nbImage = 0 nbImage = 0
if hasattr(self, 'liste_tuiles'):
self.listing_by_tile()
# print(urlencode(self.requete))
# print(self.coord_box_dd())
self.requete['box'] = self.coord_box_dd()
# sys.exit(1)
self.liste_annees = ["2019"]
self.niveau = 'LEVEL2A'
# Boucle sur les annees # Boucle sur les annees
self.logger.info("Images availables") # self.logger.info("Images availables")
for annee in self.liste_annees: for annee in self.liste_annees:
self.logger.info("=============== {0} ===============".format(annee)) self.logger.info("=============== {0} ===============".format(annee))
self.url = "{0}/{1}/api/collections/{2}/search.json?lang=fr&processingLevel={3}&_pretty=true&q={4}&box={5}&maxRecord=500".format(self.serveur, self.resto, self._capteur, self.niveau, annee, self.coord_box_dd())
if self.emprise :
# self.url = "{0}/{1}/api/collections/{2}/search.json?lang=fr&processingLevel={3}&_pretty=true&q={4}&box={5}&maxRecord=500".format(self.serveur, self.resto, self.capteur, self.niveau, annee, self.coord_box_dd())
self.url = "{0}/{1}/api/collections/{2}/search.json?{3}".format(self.serveur, self.resto, self.capteur, urlencode(self.requete))
else :
self.url = "{0}/{1}/api/collections/{2}/search.json?lang=fr&processingLevel={3}&_pretty=true&q={4}&location={5}&maxRecord=500".format(self.serveur, self.resto, self.capteur, self.niveau, annee, self.liste_tuiles[0])
# S'il existe une autre image, suivante = vrai, faux sinon # S'il existe une autre image, suivante = vrai, faux sinon
suivante = True suivante = True
print(self.url)
# Tant qu'il existe une image à traiter # Tant qu'il existe une image à traiter
while suivante: while suivante:
...@@ -212,7 +224,7 @@ class Archive(): ...@@ -212,7 +224,7 @@ class Archive():
link_archive = data_Dict['features'][d]['properties']['services']['download']['url'].replace("\\", "") link_archive = data_Dict['features'][d]['properties']['services']['download']['url'].replace("\\", "")
url_archive = link_archive.replace(self.resto, "rocket/#") url_archive = link_archive.replace(self.resto, "rocket/#")
archive_download = url_archive.replace("/download", "") # Url de l'archive à télécharger archive_download = url_archive.replace("/download", "") # Url de l'archive à télécharger
out_archive = "{0}/{1}.tgz".format(self.dossier_sortie, name_archive)# Nom de sortie de l'archive out_archive = "{0}/{1}.zip".format(self.dossier_sortie, name_archive) # Nom de sortie de l'archive
if "SENTINEL2X" not in out_archive or self.niveau == "LEVEL3A": if "SENTINEL2X" not in out_archive or self.niveau == "LEVEL3A":
self.liste_archive.append([archive_download, out_archive, feature_id]) self.liste_archive.append([archive_download, out_archive, feature_id])
...@@ -239,7 +251,7 @@ class Archive(): ...@@ -239,7 +251,7 @@ class Archive():
self.logger.info("{0} image(s) correspondent aux critères.".format(len(self.liste_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, groupe="Date"): def download_auto(self):
""" """
Méthode pour télécharger les archives sur le site Theia. Méthode pour télécharger les archives sur le site Theia.
...@@ -256,17 +268,18 @@ class Archive(): ...@@ -256,17 +268,18 @@ class Archive():
# Url pour obtenir l'authentification # Url pour obtenir l'authentification
url = "{0}/services/authenticate/".format(self.serveur) url = "{0}/services/authenticate/".format(self.serveur)
# Dictionnaire contenant les informations du proxy s'il existe # Dictionnaire contenant les informations du proxy s'il existe
proxyDict = { proxyDict = {
"http" : "{0}".format(proxy), \ "http" : "{0}".format(self.proxy), \
"https" : "{0}".format(proxy), \ "https" : "{0}".format(self.proxy), \
"ftp" : "{0}".format(proxy) \ "ftp" : "{0}".format(self.proxy) \
} }
# Dictionnaire contenant les informations d'authentification # Dictionnaire contenant les informations d'authentification
payload = { payload = {
'ident' : "{0}".format(identifiant),\ 'ident' : "{0}".format(self.id),\
'pass' : "{0}".format(mdp) 'pass' : "{0}".format(self.mdp)
} }
# Requête pour obtenir le jeton d'identification # Requête pour obtenir le jeton d'identification
...@@ -331,7 +344,6 @@ class Archive(): ...@@ -331,7 +344,6 @@ class Archive():
else : else :
sauvegarde[str(date)] = {} sauvegarde[str(date)] = {}
sauvegarde[str(date)][nom_image] = "False" sauvegarde[str(date)][nom_image] = "False"
sauvegarde[str(date)]["NDVI"] = "False"
if date in dico_date.keys() : if date in dico_date.keys() :
dico_date[date].append(l) dico_date[date].append(l)
...@@ -352,17 +364,18 @@ class Archive(): ...@@ -352,17 +364,18 @@ class Archive():
self.logger.info("Requête pour l'archive : {0}".format(img[1].split("/")[-1])) self.logger.info("Requête pour l'archive : {0}".format(img[1].split("/")[-1]))
# Url de l'archive # Url de l'archive
url = "{0}/{1}/collections/{2}/{3}/download/?issuerId=theia".format(self.serveur, self.resto, self._capteur, img[2]) url = "{0}/{1}/collections/{2}/{3}/download/?issuerId=theia".format(self.serveur, self.resto, self.capteur, img[2])
self.logger.debug("url : {}".format(url)) self.logger.debug("url : {}".format(url))
# Requête pour récupérer l'archive # Requête pour récupérer l'archive
reponse = requests.get(url, headers=head, proxies=proxyDict) reponse = requests.get(url, headers=head, proxies=proxyDict)
if extraction : if self.extraction :
# Ajout de l'archive à la liste # Ajout de l'archive à la liste
liste_content.append(reponse.content) liste_content.append(reponse.content)
del reponse del reponse
else : else :
if groupe == "Tuile" : if self.groupe == "Tuile" :
dossier_archive = "{0}/{1}/Archive".format(self.dossier_sortie, id_tuile.search(img[1]).group(0)) dossier_archive = "{0}/{1}/Archive".format(self.dossier_sortie, id_tuile.search(img[1]).group(0))
else: else:
dossier_archive = "{0}/{1}/Archive".format(self.dossier_sortie, cle[:4]) dossier_archive = "{0}/{1}/Archive".format(self.dossier_sortie, cle[:4])
...@@ -371,10 +384,11 @@ class Archive(): ...@@ -371,10 +384,11 @@ class Archive():
if not os.path.exists(dossier_archive): if not os.path.exists(dossier_archive):
os.makedirs(dossier_archive) os.makedirs(dossier_archive)
print("{0}/{1}".format(dossier_archive, img[1].split("/")[-1]))
with open("{0}/{1}".format(dossier_archive, img[1].split("/")[-1]), "wb") as fichier: with open("{0}/{1}".format(dossier_archive, img[1].split("/")[-1]), "wb") as fichier:
fichier.write(reponse.content) fichier.write(reponse.content)
if extraction : if self.extraction :
# Traitement des images (fusion, découpage selon la zone d'étude ...) # Traitement des images (fusion, découpage selon la zone d'étude ...)
self.traitement_images(cle, liste_content) self.traitement_images(cle, liste_content)
del liste_content del liste_content
...@@ -499,6 +513,7 @@ class Archive(): ...@@ -499,6 +513,7 @@ class Archive():
# Pour chaque archive # Pour chaque archive
for idx, content in enumerate(liste_content) : for idx, content in enumerate(liste_content) :
# Lecture de l'archive # Lecture de l'archive
tzip = zipfile.ZipFile(io.BytesIO(content)) tzip = zipfile.ZipFile(io.BytesIO(content))
...@@ -518,6 +533,8 @@ class Archive(): ...@@ -518,6 +533,8 @@ class Archive():
mmap_name = "/vsimem/"+uuid4().hex mmap_name = "/vsimem/"+uuid4().hex
liste_mem.append(mmap_name) liste_mem.append(mmap_name)
gdal.FileFromMemBuffer(mmap_name, tzip.read(img)) gdal.FileFromMemBuffer(mmap_name, tzip.read(img))
print(mmap_name)
sys.exit(1)
liste_bandes.append(gdal.Open(mmap_name)) liste_bandes.append(gdal.Open(mmap_name))
# On fusionne les différentes bandes en une seule image # On fusionne les différentes bandes en une seule image
......
#!/bin/env python3 #!/bin/env python3
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
import sys
import argparse import argparse
import logging import logging
from logging.handlers import RotatingFileHandler from logging.handlers import RotatingFileHandler
...@@ -9,6 +10,10 @@ from contextlib import contextmanager ...@@ -9,6 +10,10 @@ from contextlib import contextmanager
import inspect import inspect
from osgeo import gdal from osgeo import gdal
import app.Constantes as Constantes import app.Constantes as Constantes
import configparser
import ast
from datetime import datetime
import app.Satellites as Satellites
class Fusion(argparse.Action): class Fusion(argparse.Action):
def __init__(self, option_strings, dest, nargs=None, **kwargs): def __init__(self, option_strings, dest, nargs=None, **kwargs):
...@@ -119,3 +124,123 @@ def clip(image, cut, form="Mem", dst="", type_sortie=gdal.GDT_Float32, nodata="N ...@@ -119,3 +124,123 @@ def clip(image, cut, form="Mem", dst="", type_sortie=gdal.GDT_Float32, nodata="N
cutlineDSName=cut, outputType=type_sortie , format=form, dstNodata=nodata) cutlineDSName=cut, outputType=type_sortie , format=form, dstNodata=nodata)
return gdal.Warp(dst, image, options=option_clip) return gdal.Warp(dst, image, options=option_clip)
def checkDate(date):
d = date.split('-')
try:
if not date :
annee = datetime.now().year
if datetime.now().month < 10 :
mois = "0{}".format(datetime.now().month)
else :
mois = datetime.now().month
if datetime.now().day < 10 :
jour = "0{}".format(datetime.now().day)
else :
jour = datetime.now().day
elif len(d) == 1 :
annee = d[0]
mois = "01"
jour = "01"
elif len(d) == 2 :
annee = d[0]
mois = d[1]
jour = "01"
elif len(d) == 3 :
annee = d[0]
mois = d[1]
jour = d[2]
datetime(int(annee), int(mois), int(jour))
return "{annee}-{mois}-{jour}".format(annee=annee, mois=mois, jour=jour)
except ValueError:
print("Format invalide")
sys.exit(1)
except IndexError:
print("La date doit être au format 'AAAA-MM-JJ', 'AAAA-MM' ou 'AAAA'")
sys.exit(1)
def get_variable(self, config):
"""
Récupération des variables dans le fichier de configuration
"""
configfile = configparser.ConfigParser()
configfile.read(config)
self.requete = dict()
self.requete["lang"] = "fr"
self.requete["_pretty"] = "true"
self.requete["maxRecord"] = 500
# Capteur utilisé
self.capteur = configfile["satellite"]["capteur"]
self.bandes = configfile["satellite"]["bandes"]
self.requete["processingLevel"] = configfile["satellite"]["processingLevel"]
# Dossier contenant les résultats
self.dossier_sortie = configfile["sortie"]["chemin"]
try:
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.date_debut = configfile["donnees"]["date_debut"]
self.requete["startDate"] = checkDate(configfile["donnees"]["date_debut"])
except Exception as e:
raise "L'année de départ est requise."
self.requete["completionDate"] = checkDate(configfile["donnees"]["date_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 = configfile["donnees"]["chemin_emprise"]
if self.emprise :
self.zone_etude = configfile["donnees"]["chemin_zone_etude"]
if not self.zone_etude :
self.zone_etude = self.emprise
else :
try:
self.liste_tuiles = list(ast.literal_eval(configfile["donnees"]["liste_tuiles"]))
except Exception as e:
raise "Aucune emprise ni aucune tuile n'a été renseigné."
self.nombre_image = configfile["donnees"]["nombre_image"]
# Identifiant, mot de passe et proxy pour le téléchargement des images Théia
self.id = configfile["theia"]["identifiant"]
self.mdp = configfile["theia"]["mdp"]
self.proxy = configfile["theia"]["proxy"]
if self.bandes == "RGB" :
self.extent_img = Satellites.SATELLITE[self.capteur][self.requete["processingLevel"]][:-1]
else :
self.extent_img = Satellites.SATELLITE[self.capteur][self.requete["processingLevel"]]
if self.requete["processingLevel"] == "LEVEL2A" :
self.nuage = Satellites.SATELLITE[self.capteur]["NUAGE"]
else :
self.nuage = None
self.serveur = Satellites.SATELLITE[self.capteur]["serveur"]
self.resto = Satellites.SATELLITE[self.capteur]["resto"]
self.token_type = Satellites.SATELLITE[self.capteur]["token_type"]
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import os
import configparser
import datetime
import glob
from osgeo import gdal
import otbApplication as otb
from app import Archive
def str2bool(v):
return v.lower() in ("true")
class Processing(object):
def __init__(self):
pass
def i_download(self):
"""
Méthode pour télécharger les images sur le site Theia Land.
"""
if not self.annee_fin :
self.annee_fin = datetime.datetime.now().year
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, groupe=self.groupe)
self.liste_dossier = dict()
for annee in os.listdir(self.resultats):