Failed to fetch fork details. Try again later.
-
Delaigue Olivier authored7b73e330
Forked from
HYCAR-Hydro / airGR
Source project has a limited visibility.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
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
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
#!/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
class Archive():
"""
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
: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.
:type emprise: Chaîne de caractères
:param sortie: Chemin du dossier de sortie.
:type sortie: Chaîne de caractères
:param annee_debut: Année à partir de laquelle on cherche les images.
:type annee_debut: Entier
:param annee_fin: Année jusqu'à laquelle on cherche les images.
:type annee_fin: Entier
"""
def __init__(self, capteur, extensions, niveau, emprise, sortie, annee_debut, annee_fin):
"""
Créé une instance de la classe 'Archive'.
"""
self._capteur = capteur
self.niveau = niveau
self.extent_img = extensions
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.dossier_sortie = sortie
self.emprise = emprise
self.liste_archive = []
self.logger = Outils.Log("log", "Archive")
self.serveur = Satellites.SATELLITE[self._capteur]["serveur"]
self.resto = Satellites.SATELLITE[self._capteur]["resto"]
self.token_type = Satellites.SATELLITE[self._capteur]["token_type"]
self.url = ''
def utm_to_latlng(self, zone, easting, northing, northernHemisphere=True):
"""
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
:param northing: Coordonnées UTM en y.
:type northing: Flottant
: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.
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):
"""
Méthode pour obtenir les coordonnées de la zone à partir du shapefile.
:returns: Une chaîne de caractères représentant les coordonnées des coins situés \
en bas à gauche et en haut à droite.
"""
# Processus to convert the UTM shapefile
# in decimal degrees shapefile with ogr2ogr in command line
utm_outfile = "{0}/UTM_{1}".format(os.path.split(self.emprise)[0], os.path.split(self.emprise)[1])
if os.path.exists(utm_outfile):
os.remove(utm_outfile)
process_tocall = ['ogr2ogr', '-f', 'Esri Shapefile', '-t_srs', 'EPSG:32631', utm_outfile, self.emprise]
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])
return area_coord_corner
def listing(self):
"""
Méthode pour lister les images disponibles sur la plateforme 'Theia Land' correspondant à la zone.
"""
nbImage = 0
# Boucle sur les annees
self.logger.info("Images availables")
for annee in self.liste_annees:
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
# Tant qu'il existe une image à traiter
while suivante:
try :
request_headers = {"User-Agent": "Magic-browser"}
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
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
# Conversion des données sous la forme d'un dictionnaire
data_Dict = defaultdict(list)
data_Dict = UserDict(eval(new_data))
# Selection des archives à télécharger
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/#")
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
if "SENTINEL2X" not in out_archive or self.niveau == "LEVEL3A":
self.liste_archive.append([archive_download, out_archive, feature_id])
# Vérification si d'autres images sont disponibles
for link in data_Dict['properties']['links'] :
if link["title"] is 'suivant' :
self.url = link['href'].replace("\\", "")
suivante = True
break
else :
suivante = False
except Exception as e:
self.logger.error("Error connexion or error variable : {0}".format(e))
sys.exit(1)
# 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))
nbImage = 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=""):
"""
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
"""
# Url pour obtenir l'authentification
url = "{0}/services/authenticate/".format(self.serveur)
# Dictionnaire contenant les informations du proxy s'il existe
proxyDict = {
"http" : "{0}".format(proxy), \
"https" : "{0}".format(proxy), \
"ftp" : "{0}".format(proxy) \
}
# Dictionnaire contenant les informations d'authentification
payload = {
'ident' : "{0}".format(identifiant),\
'pass' : "{0}".format(mdp)
}
# Requête pour obtenir le jeton d'identification
reponse = requests.post(url, data=payload, proxies=proxyDict)
# Récupération du jeton d'identification au format texte ou json
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(.+?)_(.+?)-'
# Dictionnaire des images à traiter regrouper par date
dico_date = dict()
sauvegarde = configparser.ConfigParser()
# Pour toutes les images disponible
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)
# Récupération du nom de l'image
nom_image = l[1].split("/")[-1][:-4]
# Lecture du fichier de sauvegarde
sauvegarde.read("{}/{}/sauvegarde.ini".format(self.dossier_sortie, date[:4]))
# Si la date existe dans le fichier de sauvegarde ...
if str(date) in sauvegarde.keys() :
# ... mais pas l'image, on la rajoute dans la liste des images à traiter.
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]
# ... et que l'image exite mais n'a pas été traité, on la rajoute dans la liste des images à traiter.
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]
# Si la date n'existe pas dans le fichier de sauvegarde, on la rajoute ainsi que l'image dans la liste des images à traiter.
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]
# Sauvegarde du fichier de sauvegarde mis à jour.
with open("{}/{}/sauvegarde.ini".format(self.dossier_sortie, date[:4]), 'w') as configfile:
sauvegarde.write(configfile)
# Pour toutes les dates à traiter
for cle in sorted(dico_date):
# Liste des archives
liste_content = []
# Pour toutes les images prisent à cette date
for idx, img in enumerate(dico_date[cle]):
self.logger.info("Requête pour l'archive : {0}".format(img[1].split("/")[-1]))
# Url de l'archive
url = "{0}/{1}/collections/{2}/{3}/download/?issuerId=theia".format(self.serveur, self.resto, self._capteur, img[2])
self.logger.debug("url : {}".format(url))
# Requête pour récupérer l'archive
reponse = requests.get(url, headers=head, proxies=proxyDict)
# Ajout de l'archive à la liste
liste_content.append(reponse.content)
del reponse
# Traitement des images (fusion, découpage selon la zone d'étude ...)
self.traitement_images(cle, liste_content)
del liste_content
# Mis à jour du fichier de sauvegarde
# Lecture du fichier de sauvegarde
sauvegarde.read("{}/{}/sauvegarde.ini".format(self.dossier_sortie, cle[:4]))
# Pour toutes les images traitées à cette date
for img in dico_date[cle] :
# Image traitée passe à vrai dans le fichier de sauvegarde
sauvegarde[str(cle)][img[1].split("/")[-1][:-4]] = "True"
# Sauvegarde du fichier de sauvegarde mis à jour.
with open("{}/{}/sauvegarde.ini".format(self.dossier_sortie, cle[:4]), 'w') as configfile:
sauvegarde.write(configfile)
self.logger.info("All images have been downloaded !")
def calcul_aire(self, tuiles_image):
"""
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.
"""
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)
mask_spec = np.isin(data, [-10000, math.isnan], invert=True)
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
def ecriture_geotiff(self, dataset, chemin):
"""
Méthode pour enregistrer une image au format GTiff.
: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
"""
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) :
data = dataset.GetRasterBand(band + 1).ReadAsArray(0, 0, dataset.RasterXSize, dataset.RasterYSize).astype(np.float32)
outdata.GetRasterBand(band + 1).WriteArray(data, 0, 0)
outdata.GetRasterBand(band + 1).FlushCache()
outdata = None
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
"""
self.logger.info("Date : {0} -> {1} image(s)".format(date, len(liste_content)))
tuiles_image = []
# Options Gdal pour la fusion des bandes
options_vrt = gdal.BuildVRTOptions(separate=True)
options_translate = gdal.TranslateOptions(format="Mem")
self.logger.info("Extraction des images")
# Pour chaque archive
for idx, content in enumerate(liste_content) :
# Lecture de l'archive
tzip = zipfile.ZipFile(io.BytesIO(content))
# Liste des fichiers dans l'archive
zip_img = tzip.namelist()
# Liste contenant les images voulues
liste_bandes = []
liste_mem = []
# Pour tous les fichiers dans l'archive
for id_ext, extension in enumerate(self.extent_img) :
# Si il s'agit d'une bande voulue (R,G,B,PIR)
img = [f for f in zip_img if extension in f][0]
# On charge l'image en mémoire
mmap_name = "/vsimem/"+uuid4().hex
liste_mem.append(mmap_name)
gdal.FileFromMemBuffer(mmap_name, tzip.read(img))
liste_bandes.append(gdal.Open(mmap_name))
# On fusionne les différentes bandes en une seule image
# 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
for mmap_name in liste_mem :
gdal.Unlink(mmap_name)
liste_bandes = None
del tzip
liste_content[idx] = None
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
# Si l'aire est positive on enregistre l'image
if aire > 0.0 :
self.logger.info("Sauvegarde des images")
dossier = "{0}/{1}/Images".format(self.dossier_sortie, date[:4])
self.logger.debug("Dossier image : {0}".format(dossier))
self.ecriture_geotiff(image, "{0}/{1}.tif".format(dossier,date))
del image