diff --git a/Roseliere.py b/Roseliere.py
index ce5f1d8c65c72d1592afbdf0fd57279a4d15773e..46d5c8d67b2f24d6bbedc90341116b2b64631459 100644
--- a/Roseliere.py
+++ b/Roseliere.py
@@ -31,6 +31,7 @@
"""
import os, sys, time
+import subprocess
from _collections import defaultdict
@@ -131,6 +132,5 @@ class Roseliere(Processing):
self.logger.info("That is {0} day(s) {1}".format(nb_day_processing, time.strftime('%Hh %Mmin%S', time.gmtime(endTime - startTime))))
if __name__ == "__main__":
-
- myapp = Roseliere()
- sys.exit(myapp.run())
+ myapp = Roseliere()
+ sys.exit(myapp.run())
\ No newline at end of file
diff --git a/app/Archive.py b/app/Archive.py
index e361ec1f05803598339337f2a5305492a63e1385..51a99f37b9cad798b021ed897cf3aa4ca14b3a71 100644
--- a/app/Archive.py
+++ b/app/Archive.py
@@ -8,12 +8,12 @@
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
-#
+#
# PHYMOBAT 3.0 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
-#
+#
# You should have received a copy of the GNU General Public License
# along with PHYMOBAT 3.0. If not, see <http://www.gnu.org/licenses/>.
@@ -22,7 +22,7 @@ import math, subprocess, json
import urllib.request
import tarfile, zipfile
import requests
-import configparser
+import configparser
import datetime
from osgeo import ogr, gdal, gdal_array
@@ -34,15 +34,14 @@ from collections import defaultdict, UserDict
import app.Constantes as Constantes
import app.Satellites as Satellites
import app.Outils as Outils
-from memory_profiler import profile
class Archive():
"""
Class to list, download and unpack Theia image archive because of a shapefile (box).
This shapefile get extent of the area.
-
+
:param captor: Name of the satellite (ex: Landsat or SpotWorldHeritage ...).
-
+
Name used to the url on website Theia Land
:type captor: str
:param list_year: Processing's year (string for one year)
@@ -54,7 +53,7 @@ class Archive():
:param repertory: Name of the archive's folder
:type repertory: str
"""
-
+
def __init__(self, captor, start_year, emprise, repertory, study_area):
"""
Create a new 'Archive' instance
@@ -69,12 +68,12 @@ class Archive():
self._list_year.append(i)
self._list_year.append(now.year)
-
+
self._box = emprise
self._repertory = repertory
self.emprise = emprise
self.study_area = study_area
-
+
self.list_archive = []
self.logger = Outils.Log("log", "archive")
@@ -84,17 +83,17 @@ class Archive():
self.token_type = Satellites.SATELLITE[self._captor]["token_type"]
self.url = '' # str : complete website JSON database
-
+
self.list_img = [] # List (dim 5) to get year, month, day, path of multispectral's images and path of cloud's images
self.single_date = [] # date list without duplication
self.fichier_sauvegarde = "{0}/sauvegarde.ini".format(self._repertory)
self.sauvegarde = configparser.ConfigParser()
self.sauvegarde.read(self.fichier_sauvegarde)
-
+
def utm_to_latlng(self, zone, easting, northing, northernHemisphere=True):
"""
- Function to convert UTM to geographic coordinates
-
+ Function to convert UTM to geographic coordinates
+
@param zone: UTM zone
@type zone: int
@param easting: Coordinates UTM in x
@@ -103,60 +102,60 @@ class Archive():
@type northing: float
@param northernHemisphere: North hemisphere or not
@type northernHemisphere: boolean
-
+
:returns: tuple -- integer on the **longitude** and **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):
@@ -164,63 +163,63 @@ class Archive():
Function to get area's coordinates of shapefile
:returns: str -- **area_coord_corner** : Area coordinates corner
-
+
--> Left bottom on x, Left bottom on y, Right top on x, Right top on y
:Example:
-
+
>>> import Archive
>>> test = Archive(captor, list_year, box, folder, repertory, proxy_enabled)
>>> coor_test = test.coord_box_dd()
>>> coor_test
'45.52, 2.25, 46.71, 3.27'
"""
-
- # Processus to convert the UTM shapefile
- # in decimal degrees shapefile with ogr2ogr in command line
+
+ # 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._box)[0], os.path.split(self._box)[1])
-
+
if os.path.exists(utm_outfile):
os.remove(utm_outfile)
-
+
process_tocall = ['ogr2ogr', '-f', 'Esri Shapefile', '-t_srs', 'EPSG:32631', utm_outfile, self._box]
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)
-
+
# AdressUrl = 'http://spirit.cnes.fr/resto/Landsat/?format=html&lang=fr&q=2013&box=' + str(LL_Min[0]) + ',' + str(LL_Min[1]) + ',' + str(LL_Max[0]) + ',' + str(LL_Max[1])
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):
"""
Function to list available archive on plateform Theia Land, and on the area
"""
-
+
# Loop on the years
self.logger.info("Images availables")
for year in self._list_year:
-
+
# first_date = year.split(',')[0]
#Â Tricks to put whether a year or a date (year-month-day)
# try:
@@ -230,15 +229,15 @@ class Archive():
# except IndexError:
self.logger.info("=============== {0} ===============".format(year))
self.url = "{0}/{1}/api/collections/{2}/search.json?lang=fr&_pretty=true&q={3}&box={4}&maxRecord=500".format(self.server, self.resto, self._captor, year, self.coord_box_dd())
-
- # Initialisation variable for a next page
+
+ # Initialisation variable for a next page
# There is a next page, next = 1
# There isn't next page, next = 0
- next_ = 1
-
+ next_ = 1
+
# To know path to download images
while next_ == 1:
-
+
try :
request_headers = {"User-Agent": "Magic-browser"}
req = urllib.request.Request(self.url, headers = request_headers) # Connexion in the database
@@ -249,11 +248,11 @@ class Archive():
# Transform the data in dictionary
data_Dict = defaultdict(list)
- data_Dict = UserDict(eval(new_data))
-
+ data_Dict = UserDict(eval(new_data))
+
# Select archives to download
for d in range(len(data_Dict['features'])):
- name_archive = data_Dict['features'][d]['properties']['productIdentifier']
+ 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/#")
@@ -261,8 +260,8 @@ class Archive():
out_archive = "{0}/{1}.tgz".format(self._repertory, name_archive)# Name after download
# if len(self.list_archive) == 0 : self.list_archive.append([archive_download, out_archive, feature_id])
if "SENTINEL2X" not in out_archive :
- self.list_archive.append([archive_download, out_archive, feature_id])
-
+ self.list_archive.append([archive_download, out_archive, feature_id])
+
# Verify if there is another page (next)
for link in data_Dict['properties']['links'] :
if link["title"] is 'suivant' :
@@ -291,11 +290,11 @@ class Archive():
"https" : "{0}".format(proxy), \
"ftp" : "{0}".format(proxy) \
}
-
+
payload = {
'ident' : "{0}".format(user),\
'pass' : "{0}".format(password)
- }
+ }
reponse = requests.post(url, data=payload, proxies=proxyDict)
@@ -319,7 +318,7 @@ class Archive():
for l in self.list_archive :
date = re.search(regex, l[1]).group(2)
nom_image = l[1].split("/")[-1][:-4]
-
+
if str(date) in self.sauvegarde.keys() :
@@ -342,7 +341,7 @@ class Archive():
else :
self.sauvegarde[str(date)] = {}
self.sauvegarde[str(date)][nom_image] = "False"
-
+
if date in dico_date.keys() :
dico_date[date].append(l)
else :
@@ -351,19 +350,17 @@ class Archive():
with open(self.fichier_sauvegarde, 'w') as configfile:
self.sauvegarde.write(configfile)
-
+
for cle in sorted(dico_date):
liste_content = []
- for img in dico_date[cle]:
+ for idx, img in enumerate(dico_date[cle]):
self.logger.info("Requête pour l'archive : {0}".format(img[1].split("/")[-1]))
url = "{0}/{1}/collections/{2}/{3}/download/?issuerId=theia".format(self.server, self.resto, self._captor, img[2])
reponse = requests.get(url, headers=head, proxies=proxyDict)
-
liste_content.append(reponse.content)
del reponse
- # sys.exit()
self.pourc_cloud(cle, liste_content)
del liste_content
@@ -376,14 +373,14 @@ class Archive():
self.logger.info("All images have been downloaded !")
def calcul_aire(self, tuiles_image):
-
- option_warp = gdal.WarpOptions(format='Mem', resampleAlg="lanczos")
+
+ 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)
+ 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)
@@ -392,44 +389,40 @@ class Archive():
return mosaic_image, area
def calcul_ennuagement(self, tuiles_nuage):
-
+
option_warp = gdal.WarpOptions(format='Mem',resampleAlg="lanczos")
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).astype(np.float16)
- mask_spec = np.isin(data, [-10000, math.isnan], invert=True)
-
- # This is the same opposite False where there is 0
- mask_cloud = np.isin(data, 0)
-
- del data
+ data = mosaic_nuage.GetRasterBand(1).ReadAsArray(0, 0, cols, rows)
- # If True in cloud mask, it take spectral image else False
- cloud = np.choose(mask_cloud, (False, mask_spec))
+ nombre_pixel_zone = 0
+ nombre_pixel_null = 0
+ nombre_pixel_non_nuage = 0
- del mask_cloud
+ for donnees in data :
+ d = donnees.tolist()
+ nombre_pixel_zone += len(d)
+ nombre_pixel_null += d.count(-10000) + np.count_nonzero(np.isnan(d))
+ nombre_pixel_non_nuage += d.count(0)
- # Sum of True. True is cloud
- dist = np.sum(cloud)
+ nombre_pixel_non_null = nombre_pixel_zone - nombre_pixel_null
- del cloud
- ennuagement = 1.0 - (float(dist)/(np.sum(mask_spec)) if np.sum(mask_spec) != 0 else 0)
+ del data
- del dist
- del mask_spec
+ ennuagement = 1.0 - (float(nombre_pixel_non_nuage)/float(nombre_pixel_non_null))
# 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 mosaic_nuage, 1.0 - (float(dist)/(np.sum(mask_spec)) if np.sum(mask_spec) != 0 else 0)
return None, ennuagement
def ecriture_geotiff(self, dataset, chemin):
gdal.AllRegister()
- driver = gdal.GetDriverByName('GTiff')
+ 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())
@@ -439,7 +432,7 @@ class Archive():
dataset.GetRasterBand(band + 1).ReadAsArray(0, 0, dataset.RasterXSize, dataset.RasterYSize).astype(np.float32)/10000.0,\
0,0)
outdata.GetRasterBand(band + 1).FlushCache()
-
+
outdata = None
def pourc_cloud(self, date, liste_content):
@@ -458,8 +451,9 @@ class Archive():
for idx, content in enumerate(liste_content) :
- tzip = zipfile.ZipFile(io.BytesIO(content))
-
+ # tzip = zipfile.ZipFile(io.BytesIO(content))
+ tzip = zipfile.ZipFile(content)
+
# Images list in the archive
zip_img = tzip.namelist()
@@ -481,7 +475,7 @@ class Archive():
liste_bandes = None
- image_nuage = [f for f in zip_img if '_CLM_R1.tif' in f][0]
+ image_nuage = [f for f in zip_img if '_CLM_R1.tif' in f][0]
mmap_name = "/vsimem/"+uuid4().hex
gdal.FileFromMemBuffer(mmap_name, tzip.read(image_nuage))
@@ -495,8 +489,8 @@ class Archive():
self.logger.info("Calcul de l'aire")
image, aire = self.calcul_aire(tuiles_image)
del tuiles_image
-
- if aire > 0.0 :
+
+ if aire > 0.0 :
self.logger.info("Calcul de l'ennuagement")
nuage, ennuagement = self.calcul_ennuagement(tuiles_nuage)
self.logger.info('Ennuagement = {0}%'.format(round(ennuagement*100, 2)))
@@ -517,4 +511,4 @@ class Archive():
else :
del tuiles_nuage
- del image
\ No newline at end of file
+ del image