Archive.py 23.2 KB
Newer Older
Commandre Benjamin's avatar
Commandre Benjamin committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import os, sys, glob, re, io
import math, subprocess
import urllib.request
import zipfile
import requests
import configparser

from osgeo import ogr, gdal
import numpy as np
from uuid import uuid4

from collections import defaultdict, UserDict
import app.Satellites as Satellites
import app.Outils as Outils
Commandre Benjamin's avatar
Commandre Benjamin committed
18
import app.Constantes as Constantes
Commandre Benjamin's avatar
Commandre Benjamin committed
19
20
21

class Archive():
    """
22
23
24
25
26
        Classe pour lister, télécharger les archives via site Theia selon un shapefile représentant la zone d'étude.

        :param capteur: Nom du satellite utilisé (ex: Landsat, Sentinel2, ...).
        :type capteur: Chaîne de caractères

27
        :param niveau: Niveau de traitement voulu.
28
        :type niveau: Chaîne de caractères
29

30
31
32
33
34
        :param emprise: Chemin vers le shapefile représentant la zone d'étude.
        :type emprise: Chaîne de caractères

        :param sortie: Chemin du dossier de sortie.
        :type sortie: Chaîne de caractères
35

36
37
        :param annee_debut: Année à partir de laquelle on cherche les images.
        :type annee_debut: Entier
38

39
40
        :param annee_fin: Année jusqu'à laquelle on cherche les images.
        :type annee_fin: Entier
Commandre Benjamin's avatar
Commandre Benjamin committed
41
42
    """

Commandre Benjamin's avatar
Commandre Benjamin committed
43
    def __init__(self, capteur, bandes, niveau, emprise, zone_etude, sortie, annee_debut, annee_fin, seuil):
Commandre Benjamin's avatar
Commandre Benjamin committed
44
        """
45
            Créé une instance de la classe 'Archive'.
Commandre Benjamin's avatar
Commandre Benjamin committed
46
        """
47
48

        self._capteur = capteur
Commandre Benjamin's avatar
Commandre Benjamin committed
49
        self.niveau = niveau
50

Commandre Benjamin's avatar
Commandre Benjamin committed
51
52
53
54
55
56
57
        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"]
58
        else :
Commandre Benjamin's avatar
Commandre Benjamin committed
59
            self.nuage = None
Commandre Benjamin's avatar
Commandre Benjamin committed
60

61
        self.liste_annees = []
Commandre Benjamin's avatar
Commandre Benjamin committed
62

63
64
        for i in range(annee_debut, annee_fin) :
             self.liste_annees.append(i)
Commandre Benjamin's avatar
Commandre Benjamin committed
65

66
        self.liste_annees.append(annee_fin)
Commandre Benjamin's avatar
Commandre Benjamin committed
67

68
        self.emprise = emprise
Commandre Benjamin's avatar
Commandre Benjamin committed
69
        self.zone_etude = zone_etude
70
        self.dossier_sortie = sortie
Commandre Benjamin's avatar
Commandre Benjamin committed
71
72
        self.seuil_nuage = seuil

73
        self.liste_archive = []
Commandre Benjamin's avatar
Commandre Benjamin committed
74
75
76

        self.logger = Outils.Log("log", "Archive")

77
78
79
        self.serveur = Satellites.SATELLITE[self._capteur]["serveur"]
        self.resto = Satellites.SATELLITE[self._capteur]["resto"]
        self.token_type = Satellites.SATELLITE[self._capteur]["token_type"]
Commandre Benjamin's avatar
Commandre Benjamin committed
80
81
82
83
84

        self.url = ''

    def utm_to_latlng(self, zone, easting, northing, northernHemisphere=True):
        """
85
86
87
88
89
90
91
            Méthode pour convertir la projection UTM en coordonnées géographique

            :param zone: Zone en projection UTM.
            :type zone: Entier

            :param easting: Coordonnées UTM en x.
            :type easting: Flottant
Commandre Benjamin's avatar
Commandre Benjamin committed
92

93
94
            :param northing: Coordonnées UTM en y.
            :type northing: Flottant
Commandre Benjamin's avatar
Commandre Benjamin committed
95

96
97
98
99
            :param northernHemisphere: Vrai si la zone se situe dans l'hémisphère Nord, faux sinon.
            :type northernHemisphere: Booléen

            :returns: Tuple d'entiers représentant la longitude et la latitude.
Commandre Benjamin's avatar
Commandre Benjamin committed
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

            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))

        r0 = a * (1 - e * e) / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (3 / 2.0))
        fact1 = n0 * math.tan(phi1) / r0

        _a1 = 500000 - easting
        dd0 = _a1 / (n0 * k0)
        fact2 = dd0 * dd0 / 2

        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

        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

        return (longitude, latitude)

    def coord_box_dd(self):
        """
156
            Méthode pour obtenir les coordonnées de la zone à partir du shapefile.
Commandre Benjamin's avatar
Commandre Benjamin committed
157

158
159
            :returns: Une chaîne de caractères représentant les coordonnées des coins situés \
            en bas à gauche et en haut à droite.
Commandre Benjamin's avatar
Commandre Benjamin committed
160
161
162
163
        """

        # Processus to convert the UTM shapefile
        # in decimal degrees shapefile with ogr2ogr in command line
164
        utm_outfile = "{0}/UTM_{1}".format(os.path.split(self.emprise)[0], os.path.split(self.emprise)[1])
Commandre Benjamin's avatar
Commandre Benjamin committed
165
166
167
168

        if os.path.exists(utm_outfile):
            os.remove(utm_outfile)

169
        process_tocall = ['ogr2ogr', '-f', 'Esri Shapefile', '-t_srs', 'EPSG:32631', utm_outfile, self.emprise]
Commandre Benjamin's avatar
Commandre Benjamin committed
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
        subprocess.call(process_tocall)

        # To get shapefile extent
        # Avec import ogr
        driver = ogr.GetDriverByName('ESRI Shapefile')
        # Open shapefile
        data_source = driver.Open(utm_outfile, 0)

        if data_source is None:
            self.logger.error('Could not open file')
            sys.exit(1)

        shp_ogr = data_source.GetLayer()
        # Extent
        extent_ = shp_ogr.GetExtent() # (xMin, xMax, yMin, yMax)

        ## Close source data
        data_source.Destroy()

        # Coordinates min in decimal degrees
        LL_min = self.utm_to_latlng(31, extent_[0],  extent_[2], northernHemisphere=True)

        # Coordinates max in decimal degrees
        LL_max = self.utm_to_latlng(31, extent_[1],  extent_[3], northernHemisphere=True)

        area_coord_corner = str(LL_min[0]) + ',' + str(LL_min[1]) + ',' + str(LL_max[0]) + ',' + str(LL_max[1])
196

Commandre Benjamin's avatar
Commandre Benjamin committed
197
198
199
200
        return area_coord_corner

    def listing(self):
        """
201
            Méthode pour lister les images disponibles sur la plateforme 'Theia Land' correspondant à la zone.
Commandre Benjamin's avatar
Commandre Benjamin committed
202
203
204
        """

        nbImage = 0
205
206

        # Boucle sur les annees
Commandre Benjamin's avatar
Commandre Benjamin committed
207
        self.logger.info("Images availables")
208
        for annee in self.liste_annees:
Commandre Benjamin's avatar
Commandre Benjamin committed
209

210
211
212
213
            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())
            # S'il existe une autre image, suivante = vrai, faux sinon
            suivante = True
Commandre Benjamin's avatar
Commandre Benjamin committed
214

215
216
            # Tant qu'il existe une image à traiter
            while suivante:
Commandre Benjamin's avatar
Commandre Benjamin committed
217
218
219

                try :
                    request_headers = {"User-Agent": "Magic-browser"}
220
221
                    req = urllib.request.Request(self.url, headers = request_headers) # Connexion à la base de données
                    data = urllib.request.urlopen(req).read() # Lecture de la base de données
Commandre Benjamin's avatar
Commandre Benjamin committed
222

223
224
                    new_data = re.sub(b"null", b"'null'", data) # Suppression de l'expression "null" pouvant causer une erreur Python
                    new_data = re.sub(b"false", b"False", new_data) # Renomme le booléen selon la syntaxe Python
Commandre Benjamin's avatar
Commandre Benjamin committed
225

226
                    # Conversion des données sous la forme d'un dictionnaire
Commandre Benjamin's avatar
Commandre Benjamin committed
227
228
229
                    data_Dict = defaultdict(list)
                    data_Dict = UserDict(eval(new_data))

230
                    # Selection des archives à télécharger
Commandre Benjamin's avatar
Commandre Benjamin committed
231
232
233
234
235
                    for d in range(len(data_Dict['features'])):
                        name_archive = data_Dict['features'][d]['properties']['productIdentifier']
                        feature_id = data_Dict["features"][d]['id']
                        link_archive = data_Dict['features'][d]['properties']['services']['download']['url'].replace("\\", "")
                        url_archive = link_archive.replace(self.resto, "rocket/#")
236
237
238
                        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

Commandre Benjamin's avatar
Commandre Benjamin committed
239
                        if "SENTINEL2X" not in out_archive or self.niveau == "LEVEL3A":
240
                            self.liste_archive.append([archive_download, out_archive, feature_id])
Commandre Benjamin's avatar
Commandre Benjamin committed
241

242
                    # Vérification si d'autres images sont disponibles
Commandre Benjamin's avatar
Commandre Benjamin committed
243
244
245
                    for link in data_Dict['properties']['links'] :
                        if link["title"] is 'suivant' :
                            self.url = link['href'].replace("\\", "")
246
                            suivante = True
Commandre Benjamin's avatar
Commandre Benjamin committed
247
248
                            break
                        else :
249
                            suivante = False
Commandre Benjamin's avatar
Commandre Benjamin committed
250
251
252
253
254

                except Exception as e:
                    self.logger.error("Error connexion or error variable : {0}".format(e))
                    sys.exit(1)

255
256
257
258
            # Si aucune image n'est disponible pour une année, le dossier correspondant n'est pas créé
            if nbImage != len(self.liste_archive) :
                if not os.path.exists("{}/{}/Images".format(self.dossier_sortie, annee)) :
                    os.makedirs("{}/{}/Images".format(self.dossier_sortie, annee))
Commandre Benjamin's avatar
Commandre Benjamin committed
259

260
            nbImage = len(self.liste_archive)
Commandre Benjamin's avatar
Commandre Benjamin committed
261

262
        self.logger.info("{0} image(s) correspondent aux critères.".format(len(self.liste_archive)))
Commandre Benjamin's avatar
Commandre Benjamin committed
263

264
    def download_auto(self, identifiant, mdp, proxy="", extraction=True):
Commandre Benjamin's avatar
Commandre Benjamin committed
265
        """
266
267
268
269
270
271
272
273
274
275
            Méthode pour télécharger les archives sur le site Theia.

            :param identifiant: Identifant pour le site Theia.
            :type identifiant: Chaîne de caractères

            :param mdp: Mot de passe pour le site Theia.
            :type mdp: Chaîne de caractères

            :param proxy: Information du proxy.
            :type proxy: Chaîne de caractères
Commandre Benjamin's avatar
Commandre Benjamin committed
276
277
        """

278
279
280
281
        # Url pour obtenir l'authentification
        url = "{0}/services/authenticate/".format(self.serveur)

        # Dictionnaire contenant les informations du proxy s'il existe
Commandre Benjamin's avatar
Commandre Benjamin committed
282
283
284
285
286
287
        proxyDict = {
            "http"  : "{0}".format(proxy), \
            "https" : "{0}".format(proxy), \
            "ftp"   : "{0}".format(proxy) \
        }

288
        # Dictionnaire contenant les informations d'authentification
Commandre Benjamin's avatar
Commandre Benjamin committed
289
        payload = {
290
291
            'ident' : "{0}".format(identifiant),\
            'pass' : "{0}".format(mdp)
Commandre Benjamin's avatar
Commandre Benjamin committed
292
293
        }

294
        # Requête pour obtenir le jeton d'identification
Commandre Benjamin's avatar
Commandre Benjamin committed
295
296
        reponse = requests.post(url, data=payload, proxies=proxyDict)

297
        # Récupération du jeton d'identification au format texte ou json
Commandre Benjamin's avatar
Commandre Benjamin committed
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
        try:
            if "text" in reponse.headers["Content-Type"] :
                token = reponse.text
            elif "json" in reponse.headers["Content-Type"] :
                token = reponse.json()["access_token"]
            else :
                raise(BaseException('Format non traité.'))
        except Exception as e:
            self.logger.error("Error during the identification request")
            sys.exit(-1)

        head = {"Authorization": "Bearer {0}".format(token)}

        regex = 'SENTINEL(.+?)_(.+?)-'

313
        # Dictionnaire des images à traiter regrouper par date
Commandre Benjamin's avatar
Commandre Benjamin committed
314
315
316
317
        dico_date = dict()

        sauvegarde = configparser.ConfigParser()

318
319
320
321
        # Pour toutes les images disponible
        for l in self.liste_archive :

            # Récupération de la date quand l'image a été prise
Commandre Benjamin's avatar
Commandre Benjamin committed
322
            date = re.search(regex, l[1]).group(2)
323
324

            # Récupération du nom de l'image
Commandre Benjamin's avatar
Commandre Benjamin committed
325
326
            nom_image = l[1].split("/")[-1][:-4]

327
328
            # Lecture du fichier de sauvegarde
            sauvegarde.read("{}/{}/sauvegarde.ini".format(self.dossier_sortie, date[:4]))
Commandre Benjamin's avatar
Commandre Benjamin committed
329

330
            # Si la date existe dans le fichier de sauvegarde ...
Commandre Benjamin's avatar
Commandre Benjamin committed
331
332
            if str(date) in sauvegarde.keys() :

333
                # ... mais pas l'image, on la rajoute dans la liste des images à traiter.
Commandre Benjamin's avatar
Commandre Benjamin committed
334
335
336
337
338
339
340
341
342
343
                if nom_image not in sauvegarde[str(date)].keys() :

                    sauvegarde[str(date)][nom_image] = "False"
                    sauvegarde[str(date)]["NDVI"] = "False"

                    if date in dico_date.keys() :
                        dico_date[date].append(l)
                    else :
                        dico_date[date] = [l]

344
                # ... et que l'image exite mais n'a pas été traité, on la rajoute dans la liste des images à traiter.
Commandre Benjamin's avatar
Commandre Benjamin committed
345
346
347
348
349
350
351
352
353
                elif sauvegarde[str(date)][nom_image] == "False" :

                    sauvegarde[str(date)]["NDVI"] = "False"

                    if date in dico_date.keys() :
                        dico_date[date].append(l)
                    else :
                        dico_date[date] = [l]

354
            # Si la date n'existe pas dans le fichier de sauvegarde, on la rajoute ainsi que l'image dans la liste des images à traiter.
Commandre Benjamin's avatar
Commandre Benjamin committed
355
356
357
358
359
360
361
362
363
364
            else :
                sauvegarde[str(date)] = {}
                sauvegarde[str(date)][nom_image] = "False"
                sauvegarde[str(date)]["NDVI"] = "False"

                if date in dico_date.keys() :
                    dico_date[date].append(l)
                else :
                    dico_date[date] = [l]

365
366
            # Sauvegarde du fichier de sauvegarde mis à jour.
            with open("{}/{}/sauvegarde.ini".format(self.dossier_sortie, date[:4]), 'w') as configfile:
Commandre Benjamin's avatar
Commandre Benjamin committed
367
368
                sauvegarde.write(configfile)

369
        # Pour toutes les dates à traiter
Commandre Benjamin's avatar
Commandre Benjamin committed
370
        for cle in sorted(dico_date):
371
            # Liste des archives
Commandre Benjamin's avatar
Commandre Benjamin committed
372
            liste_content = []
373
374

            # Pour toutes les images prisent à cette date
Commandre Benjamin's avatar
Commandre Benjamin committed
375
            for idx, img in enumerate(dico_date[cle]):
376

Commandre Benjamin's avatar
Commandre Benjamin committed
377
                self.logger.info("Requête pour l'archive : {0}".format(img[1].split("/")[-1]))
378
379
                # Url de l'archive
                url = "{0}/{1}/collections/{2}/{3}/download/?issuerId=theia".format(self.serveur, self.resto, self._capteur, img[2])
Commandre Benjamin's avatar
Commandre Benjamin committed
380
                self.logger.debug("url : {}".format(url))
381
                # Requête pour récupérer l'archive
Commandre Benjamin's avatar
Commandre Benjamin committed
382
383
                reponse = requests.get(url, headers=head, proxies=proxyDict)

384
385
386
387
388
389
390
391
392
393
394
395
396
                if extraction :
                    # Ajout de l'archive à la liste
                    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:
                        fichier.write(reponse.content)

            if extraction :
                # Traitement des images (fusion, découpage selon la zone d'étude ...)
                self.traitement_images(cle, liste_content)
                del liste_content
Commandre Benjamin's avatar
Commandre Benjamin committed
397

398
399
400
            # Mis à jour du fichier de sauvegarde
            # Lecture du fichier de sauvegarde
            sauvegarde.read("{}/{}/sauvegarde.ini".format(self.dossier_sortie, cle[:4]))
Commandre Benjamin's avatar
Commandre Benjamin committed
401

402
            # Pour toutes les images traitées à cette date
Commandre Benjamin's avatar
Commandre Benjamin committed
403
            for img in dico_date[cle] :
404
                # Image traitée passe à vrai dans le fichier de sauvegarde
Commandre Benjamin's avatar
Commandre Benjamin committed
405
406
                sauvegarde[str(cle)][img[1].split("/")[-1][:-4]] = "True"

407
408
            # Sauvegarde du fichier de sauvegarde mis à jour.
            with open("{}/{}/sauvegarde.ini".format(self.dossier_sortie, cle[:4]), 'w') as configfile:
Commandre Benjamin's avatar
Commandre Benjamin committed
409
410
411
412
                sauvegarde.write(configfile)

        self.logger.info("All images have been downloaded !")

Commandre Benjamin's avatar
Commandre Benjamin committed
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
    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)))
444

Commandre Benjamin's avatar
Commandre Benjamin committed
445
446
447
448
449
        # 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


Commandre Benjamin's avatar
Commandre Benjamin committed
450
    def calcul_aire(self, tuiles_image):
451
452
453
454
455
456
457
458
        """
            Méthode pour fusionner un ensemble d'images et calculer l'aire de l'image obtenue.

            :param tuiles_image: Liste des images à fusionner.
            :type tuiles_image: Liste d'images

            :returns: Un tuple contenant l'image fusionnée ainsi que son aire.
        """
Commandre Benjamin's avatar
Commandre Benjamin committed
459
460
461
462
463
464
465
466

        option_warp     = gdal.WarpOptions(format='Mem', resampleAlg="lanczos")
        mosaic_image    = gdal.Warp("", tuiles_image, options=option_warp)

        rows = mosaic_image.RasterYSize # Rows number
        cols = mosaic_image.RasterXSize # Columns number

        data        = mosaic_image.GetRasterBand(1).ReadAsArray(0, 0, cols, rows).astype(np.float16)
Commandre Benjamin's avatar
Commandre Benjamin committed
467
        mask_spec   = ~(np.isnan(data) | np.isin(data, [-10000]))
Commandre Benjamin's avatar
Commandre Benjamin committed
468
469

        area = float((np.sum(mask_spec) * mosaic_image.GetGeoTransform()[1] * abs(mosaic_image.GetGeoTransform()[-1]) )/10000)
Commandre Benjamin's avatar
Commandre Benjamin committed
470
        # self.logger.info("Aire = {0} ha".format(area))
Commandre Benjamin's avatar
Commandre Benjamin committed
471
472
473
474

        return mosaic_image, area

    def ecriture_geotiff(self, dataset, chemin):
475
476
        """
            Méthode pour enregistrer une image au format GTiff.
Commandre Benjamin's avatar
Commandre Benjamin committed
477

478
479
480
481
482
483
            :param dataset: Données de l'image.
            :type dataset: GDALDataset

            :param chemin: Chemin de l'image de sortie.
            :type chemin: Chaîne de caractères
        """
Commandre Benjamin's avatar
Commandre Benjamin committed
484
485
486
487
488
489
490
491
        gdal.AllRegister()

        driver = gdal.GetDriverByName('GTiff')
        outdata = driver.Create(chemin, dataset.RasterXSize, dataset.RasterYSize, dataset.RasterCount, gdal.GDT_Float32)
        outdata.SetGeoTransform(dataset.GetGeoTransform())
        outdata.SetProjection(dataset.GetProjection())

        for band in range(dataset.RasterCount) :
492
493
            data = dataset.GetRasterBand(band + 1).ReadAsArray(0, 0, dataset.RasterXSize, dataset.RasterYSize).astype(np.float32)
            outdata.GetRasterBand(band + 1).WriteArray(data, 0, 0)
Commandre Benjamin's avatar
Commandre Benjamin committed
494
495
496
497
            outdata.GetRasterBand(band + 1).FlushCache()

        outdata = None

498
499
500
501
502
503
504
505
506
507
508
509
510
511
    def traitement_images(self, date, liste_content):
        """
            Pour une date donnée, extrait pour chaque archive correspondant à cette date
            les images correspondant aux bandes rouge, verte, bleue et proche infrarouge,
            les combines en une seule image, effectue la mosaïque des images,
            coupe selon l'emprise de la zone d'étude puis, si l'aire de l'image est supérieure à 0 l'enregistre.

            :param date: Date à laquelle ont été prises les images.
            :type date: Entier

            :param liste_content: Liste des archives contenant les images.
            :type liste_content: Liste d'octets

        """
Commandre Benjamin's avatar
Commandre Benjamin committed
512
513
514
515

        self.logger.info("Date : {0} -> {1} image(s)".format(date, len(liste_content)))

        tuiles_image = []
Commandre Benjamin's avatar
Commandre Benjamin committed
516
        tuiles_nuage = []
Commandre Benjamin's avatar
Commandre Benjamin committed
517

518
        # Options Gdal pour la fusion des bandes
519
        options_vrt = gdal.ParseCommandLine('-resolution highest -srcnodata -10000 -vrtnodata NaN -separate')
Commandre Benjamin's avatar
Commandre Benjamin committed
520
521
522

        self.logger.info("Extraction des images")

523
        # Pour chaque archive
Commandre Benjamin's avatar
Commandre Benjamin committed
524
        for idx, content in enumerate(liste_content) :
525
            # Lecture de l'archive
Commandre Benjamin's avatar
Commandre Benjamin committed
526
527
            tzip = zipfile.ZipFile(io.BytesIO(content))

528
            # Liste des fichiers dans l'archive
Commandre Benjamin's avatar
Commandre Benjamin committed
529
530
            zip_img = tzip.namelist()

531
            # Liste contenant les images voulues
Commandre Benjamin's avatar
Commandre Benjamin committed
532
            liste_bandes = []
533

Commandre Benjamin's avatar
Commandre Benjamin committed
534
            liste_mem = []
535
            # Pour tous les fichiers dans l'archive
536
            for id_ext, extension in enumerate(self.extent_img) :
537
                # Si il s'agit d'une bande voulue (R,G,B,PIR)
Commandre Benjamin's avatar
Commandre Benjamin committed
538
539
                img = [f for f in zip_img if extension in f][0]

540
                # On charge l'image en mémoire
Commandre Benjamin's avatar
Commandre Benjamin committed
541
542
543
544
545
                mmap_name = "/vsimem/"+uuid4().hex
                liste_mem.append(mmap_name)
                gdal.FileFromMemBuffer(mmap_name, tzip.read(img))
                liste_bandes.append(gdal.Open(mmap_name))

546
547
548
            # On fusionne les différentes bandes en une seule image
            # on découpe l'image selon l'emprise
            # et on conserve l'image obtenue
Commandre Benjamin's avatar
Commandre Benjamin committed
549
550
            vrt = gdal.BuildVRT("", liste_bandes, options=options_vrt)

Commandre Benjamin's avatar
Commandre Benjamin committed
551
552
553
            image = Outils.clip(vrt, self.zone_etude)
            _, aire = self.calcul_aire(image)

554
555
            del image

Commandre Benjamin's avatar
Commandre Benjamin committed
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
            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))

Commandre Benjamin's avatar
Commandre Benjamin committed
573
574
575
            for mmap_name in liste_mem :
                gdal.Unlink(mmap_name)

Commandre Benjamin's avatar
Commandre Benjamin committed
576
577
            liste_mem = []

Commandre Benjamin's avatar
Commandre Benjamin committed
578
579
580
581
582
583
584
            liste_bandes = None

            del tzip
            liste_content[idx] = None

        del liste_content

585
        if len(tuiles_image) > 0 :
Commandre Benjamin's avatar
Commandre Benjamin committed
586
587
588
589

            if len(tuiles_nuage) == 0 or self.calcul_ennuagement(tuiles_nuage) <= self.seuil_nuage:

                self.logger.info("Sauvegarde des images")
Commandre Benjamin's avatar
Commandre Benjamin committed
590

Commandre Benjamin's avatar
Commandre Benjamin committed
591
                dossier = "{0}/{1}/Images".format(self.dossier_sortie, date[:4])
Commandre Benjamin's avatar
Commandre Benjamin committed
592

Commandre Benjamin's avatar
Commandre Benjamin committed
593
                self.logger.debug("Dossier image : {0}".format(dossier))
Commandre Benjamin's avatar
Commandre Benjamin committed
594

595
596
                # On effectue une mosaïque des images qu'on découpe selon l'emprise.
                Outils.clip(tuiles_image, self.emprise, form="GTiff", dst="{0}/{1}.tif".format(dossier,date))
Commandre Benjamin's avatar
Commandre Benjamin committed
597

Commandre Benjamin's avatar
Commandre Benjamin committed
598
                del tuiles_nuage
Commandre Benjamin's avatar
Commandre Benjamin committed
599

600
            del tuiles_image