diff --git a/PHYMOBAT.py b/PHYMOBAT.py
index 3a4db08e17f6e22d364b7e263d02ca329198f5dc..972037145e9c508d8547395aa44bb04e84060495 100644
--- a/PHYMOBAT.py
+++ b/PHYMOBAT.py
@@ -42,7 +42,6 @@ except :
import webbrowser
import lxml.etree as ET
-import os.path
import Popup
from Processing import Processing
@@ -81,8 +80,8 @@ class PHYMOBAT(QtWidgets.QMainWindow, Processing):
# Chargement et lancement automatique pour les tests à supprimer
# self.open_backup("/media/commandre/8ce2dd88-1425-4161-a3c1-47f566c0285a/data_test/2018.xml")
- self.open_backup("/media/commandre/8ce2dd88-1425-4161-a3c1-47f566c0285a/Gironde_Test/gironde_total.xml")
- self.ok_button()
+ # self.open_backup("/media/commandre/8ce2dd88-1425-4161-a3c1-47f566c0285a/Gironde_Test/gironde_total.xml")
+ # self.ok_button()
def initUI(self):
@@ -881,6 +880,7 @@ class PHYMOBAT(QtWidgets.QMainWindow, Processing):
# Another check box to launch VHRS texture processing. If checked, vs = 1.
if self.ui.checkBox_VHRS.isChecked():
vs = 1
+
self.i_images_processing(vs) #Â function to launch the image processing
# To launch texture processing only
@@ -934,7 +934,7 @@ class PHYMOBAT(QtWidgets.QMainWindow, Processing):
# Compute a output slope raster
if self.path_mnt :
self.i_slope()
-
+
# Save the sample features
self.add_sample()
@@ -975,12 +975,12 @@ class PHYMOBAT(QtWidgets.QMainWindow, Processing):
# in_backup = QtWidgets.QFileDialog.getOpenFileName(self, "Open backup", os.getcwd(), '*.xml')[0]
- if test is None :
+ if test is False :
in_backup = QtWidgets.QFileDialog.getOpenFileName(self, "Open backup", os.getcwd(), '*.xml')[0]
else :
in_backup = test
- if in_backup != None and len(in_backup) > 0 :
+ if in_backup and len(in_backup) > 0 :
# Parse the xml file if the user choose a file
tree = ET.parse(str(in_backup))
diff --git a/Processing.py b/Processing.py
index 7d10993bc9bb55bdfc6153e841efc0c2dc6570be..96a26699a773222b7a5a621185f182dcc83c50f7 100644
--- a/Processing.py
+++ b/Processing.py
@@ -17,7 +17,7 @@
# You should have received a copy of the GNU General Public License
# along with PHYMOBAT 2.0. If not, see <http://www.gnu.org/licenses/>.
-import os, sys, math
+import os, sys, math, psutil, glob
import numpy as np
import subprocess
from sklearn.ensemble import RandomForestClassifier
@@ -281,8 +281,6 @@ class Processing(object):
check_L8.mosaic_by_date(date)
# Search cloud's percentage, select image and compute ndvi index
- # if cloud_pourcentage < CLOUD_THRESHOLD
-
if check_L8.pourc_cloud() < Constantes.CLOUD_THRESHOLD :
tmp = date
check_L8.calcul_ndvi()
@@ -292,25 +290,36 @@ class Processing(object):
# Compute temporal stats on ndvi index [min, max, std, min-max]
spectral_trans = np.transpose(np.array(spectral_out, dtype=object))
- stats_name = ['Min', 'Date', 'Max', 'Std', 'MaxMin']
- stats_ndvi, stats_cloud = current_list.calc_serie_stats(spectral_trans)
+ del spectral_out
+
+ # stats_name = ['Min', 'Date', 'Max', 'Std', 'MaxMin']
+ stats_name = ['Min', 'Max']
+ path_stat = current_list.calc_serie_stats(spectral_trans, stats_name, self.folder_processing)
+
+ for idx, i in enumerate(stats_name):
+ self.out_ndvistats_folder_tab[i] = path_stat[idx]
- # Stats cloud raster
- stats_cloud_raster = "{0}/Cloud_number_{1}.TIF".format(self.folder_processing, self.captor_project)
+ # self.logger.debug(self.out_ndvistats_folder_tab)
+ # raise(BaseException("Tst"))
- check_L8.complete_raster(\
- *check_L8.create_empty_raster(stats_cloud_raster, stats_cloud,\
- check_L8.raster_data(self.check_download.list_img[0][4])[1]),\
- stats_cloud)
+ # stats_ndvi, stats_cloud = current_list.calc_serie_stats(spectral_trans, self.folder_processing)
- # Stats ndvi rasters
- for stats_index in range(len(stats_ndvi)):
- out_ndvistats_raster = "{0}/{1}_{2}.TIF".format(self.folder_processing, stats_name[stats_index], self.captor_project)
- self.out_ndvistats_folder_tab[stats_index] = out_ndvistats_raster
- check_L8.complete_raster(\
- *check_L8.create_empty_raster(out_ndvistats_raster, stats_ndvi[stats_index],\
- check_L8.raster_data(self.check_download.list_img[0][4])[1]),\
- stats_ndvi[stats_index])
+ # # Stats cloud raster
+ # stats_cloud_raster = "{0}/Cloud_number_{1}.TIF".format(self.folder_processing, self.captor_project)
+
+ # check_L8.complete_raster(\
+ # *check_L8.create_empty_raster(stats_cloud_raster, stats_cloud,\
+ # check_L8.raster_data(self.check_download.list_img[0][4])[1]),\
+ # stats_cloud)
+
+ # # Stats ndvi rasters
+ # for stats_index in range(len(stats_ndvi)):
+ # out_ndvistats_raster = "{0}/{1}_{2}.TIF".format(self.folder_processing, stats_name[stats_index], self.captor_project)
+ # self.out_ndvistats_folder_tab[stats_index] = out_ndvistats_raster
+ # check_L8.complete_raster(\
+ # *check_L8.create_empty_raster(out_ndvistats_raster, stats_ndvi[stats_index],\
+ # check_L8.raster_data(self.check_download.list_img[0][4])[1]),\
+ # stats_ndvi[stats_index])
check_L8.logger.close()
current_list.logger.close()
@@ -350,11 +359,12 @@ class Processing(object):
current_path_ortho.image = self.path_ortho
current_path_ortho.vect = self.path_area
current_path_ortho.output = "{0}/MosaiqueOrtho/Clip_{1}".format(self.folder_processing, os.path.basename(self.path_ortho))
-
+
+ self.path_ortho = current_path_ortho.output
path_ortho = current_path_ortho.clip_raster()
texture_irc = Vhrs(path_ortho, self.mp)
- self.out_ndvistats_folder_tab['sfs'] = texture_irc.out_sfs
+ self.out_ndvistats_folder_tab['sfs'] = texture_irc.out_sfs
self.out_ndvistats_folder_tab['haralick'] = texture_irc.out_haralick
current_path_ortho.logger.close()
@@ -369,37 +379,34 @@ class Processing(object):
######################### Multiprocessing ? #########################
+ # self.i_vhrs()
+ # sys.exit(0)
+
mgr = BaseManager()
mgr.register('defaultdict', defaultdict, DictProxy)
mgr.start()
- self.out_ndvistats_folder_tab = mgr.defaultdict(list)
-
- self.i_vhrs()
- self.i_img_sat()
+ self.out_ndvistats_folder_tab = mgr.defaultdict(list)
-
- # p_img_sat = Process(target=self.i_img_sat)
- # p_vhrs = Process(target=self.i_vhrs)
+ p_img_sat = Process(target=self.i_img_sat)
+ p_vhrs = Process(target=self.i_vhrs)
- # self.logger.debug("Multiprocessing : {0}".format(self.mp))
+ self.logger.debug("Multiprocessing : {0}".format(self.mp))
- # if self.mp == Constantes.MULTIPROCESSING_DISABLE:
- # p_img_sat.start()
- # p_img_sat.join()
- # p_vhrs.start()
- # p_vhrs.join()
- # else :
- # p_img_sat.start()
- # p_vhrs.start()
- # p_img_sat.join()
- # p_vhrs.join()
-
- # raise('i_images_processing')
+ if self.mp == Constantes.MULTIPROCESSING_DISABLE:
+ p_img_sat.start()
+ p_img_sat.join()
+ p_vhrs.start()
+ p_vhrs.join()
+ else :
+ p_img_sat.start()
+ p_vhrs.start()
+ p_img_sat.join()
+ p_vhrs.join()
###################################################################
- self.raster_path["Min"]["PATH"] = self.out_ndvistats_folder_tab[0]
+ self.raster_path["Min"]["PATH"] = self.out_ndvistats_folder_tab['Min']
self.raster_path["Min"]["BAND"] = 1
for i in range(1,7) :
@@ -415,7 +422,7 @@ class Processing(object):
self.raster_path["MNT"]["PATH"] = self.path_mnt
self.raster_path["MNT"]["BAND"] = 1
- self.raster_path["MAX"]["PATH"] = self.out_ndvistats_folder_tab[2]
+ self.raster_path["MAX"]["PATH"] = self.out_ndvistats_folder_tab['Max']
self.raster_path["MAX"]["BAND"] = 1
self.logger.info("End of images processing !")
@@ -566,10 +573,10 @@ class Processing(object):
# Add the validation shapefile
self.valid_shp.append([sample_rd[sple].vector_val, kwargs['fieldname'], kwargs['class']])
- self.logger.debug(self.valid_shp)
-
items = list(self.raster_path.items())
+ self.logger.debug("Raster path : {0}".format(self.raster_path.items()))
+
for lbo in range(len(self.raster_path)):
self.logger.debug("lbo : {0}".format(lbo))
kwargs['rank'] = lbo
@@ -581,6 +588,9 @@ class Processing(object):
# To convert the dictionnary in a list
# for key, value in sample_rd[sple].stats_dict.iteritems():
for key, value in sample_rd[sple].stats_dict.items():
+
+ print("{0} : {1}".format(key, value))
+
# X_rf.append([-10000 if (math.isnan(x) or math.isinf(x)) else x for x in value])
# To set the grassland class of the RPG and PIAO (same class)
if sple == 2 :
@@ -662,15 +672,19 @@ class Processing(object):
# Stats backup file
file_stats = "{0}/Stats_raster_spectral_texture.stats".format(os.path.dirname(items[ind_th][1]["PATH"]))
- if not os.path.exists(file_stats):
- exist_stats = 0 # The sats file doesn't exist
+ """
+ TODO : Regarder si le fichier stats contient les infos pas seulement
+ tester l'existence
+ """
+ if not os.path.exists(file_stats) : # or True:
+ exist_stats = 0# The sats file doesn't exist
# Open a stats backup to avoid computing again (Gain of time)
f_out = open(file_stats, "w")
p = []
kwargs = {}
X_out_rf = [] # Variable list to compute decision tree with random forest method
- if exist_stats == 0:
+ if exist_stats == 0 :
out_carto.stats_dict = mgr.defaultdict(list)
for i in range(len(multi_process_var)) :
@@ -704,7 +718,7 @@ class Processing(object):
index_in_stats = index_in_stats + 1
out_carto.stats_dict[index_in_stats] = eval(x_in.strip('\n'))
X_out_rf.append(eval(x_in.strip('\n')))
-
+
predicted_rf = self.rf.predict(X_out_rf)
self.logger.debug("Taille sample_name : {0}".format(len(self.sample_name)))
@@ -756,7 +770,10 @@ class Processing(object):
rpg_tif.clip_vector(os.path.dirname(self.sample_name[0]))
rpg_tif.vector_data()
+ # raise(BaseException("T2"))
+
kwargs['choice_nb_b'] = 1
+ self.logger.debug(self.path_ortho)
out_carto.stats_rpg_tif = out_carto.zonal_stats_pp(rpg_tif.layer_rasterization(self.path_ortho, 'CODE_GROUP', **kwargs))
# Final cartography
@@ -920,7 +937,7 @@ class Processing(object):
opt['class'] = str(class_validate) # Add integer classes
# Set the new shapefile
- val[0] = val[0][:-4] + '_.shp'
+ val[0] = "{0}_.shp".format(val[0][:-4])
val[1] = opt['fieldname']
val[2] = opt['class']
@@ -951,7 +968,6 @@ class Processing(object):
valid.ex_raster = list(self.raster_path.items())[0][1]["PATH"]
# TODO: Call the RasterSat_by_Date class here instead of the Precision_moba class
-
self.logger.debug("Name moba : {0}".format(self.output_name_moba))
valid.preprocess_to_raster_precision(self.output_name_moba, 'FBPHY_CODE') # To the classification's data
diff --git a/RasterSat_by_date.py b/RasterSat_by_date.py
index 84ee4ab2f4a2601030d94b34b614e4c5dc02eb13..d87702ca8100c4e0c9eb5d54a3a1632fa4036e73 100644
--- a/RasterSat_by_date.py
+++ b/RasterSat_by_date.py
@@ -21,6 +21,7 @@ import os, sys
import math, subprocess
import osgeo.gdal
+import otbApplication as otb
import Outils
import Satellites
@@ -179,7 +180,8 @@ class RasterSat_by_date(object):
def raster_data(self, img):
"""
Function to extract raster information.
- Return table of pixel values and raster information like line number, pixel size, ... (gdal pointer)
+ Return table of pixel values and raster information like
+ line number, pixel size, ... (gdal pointer)
@param img : Raster path
@type img : str
@@ -192,7 +194,7 @@ class RasterSat_by_date(object):
# Load Gdal's drivers
osgeo.gdal.AllRegister()
- self.logger.debug("Img : {0}".format(img))
+ self.logger.debug("Image : {0}".format(img))
# Loading input raster
self.logger.info('Loading input raster : {0}'.format(os.path.basename(img)))
@@ -214,16 +216,19 @@ class RasterSat_by_date(object):
# Table's declaration
data = []
+
for band in range(nbband):
canal = in_ds.GetRasterBand(band + 1) # Select a band
if nbband == 1:
- data = canal.ReadAsArray(0, 0, cols, rows).astype(np.float32) # Assign pixel values at th0e data
+ # Assign pixel values at the data
+ data = canal.ReadAsArray(0, 0, cols, rows).astype(np.float16)
else:
- data.append(canal.ReadAsArray(0, 0, cols, rows).astype(np.float32))
+ data.append(canal.ReadAsArray(0, 0, cols, rows).astype(np.float16))
###################################
+
# Close input raster
_in_ds = in_ds
in_ds = None
@@ -250,7 +255,7 @@ class RasterSat_by_date(object):
data_cloud, info_cloud = self.raster_data(self.cloudiness_pourcentage[1])
# Add cloud pixel value table
- self.cloudiness_pourcentage.append(data_cloud)
+ # self.cloudiness_pourcentage.append(data_cloud)
# Extent of the images :
# ex : array([ True, True, True, True, False, True, True, True, True], dtype=bool)
@@ -273,44 +278,39 @@ class RasterSat_by_date(object):
dist = np.sum(cloud)
# Computer cloud's percentage with dist (sum of cloud) by sum of the image's extent
- try :
- nb0 = float(dist)/(np.sum(mask_spec))
- self.logger.info('For {0}, cloudy cover = {1}%'.format(os.path.basename(self.cloudiness_pourcentage[0]), 100 - round(nb0*100, 2) ) )
- except ZeroDivisionError:
- nb0 = 0
- self.logger.info("The raster isn\'t in the area !")
+ nb0 = float(dist)/(np.sum(mask_spec)) if np.sum(mask_spec) != 0 else 0
+ self.logger.info('For {0}, cloudy cover = {1}%'.format(os.path.basename(self.cloudiness_pourcentage[0]), 100 - round(nb0*100, 2)))
- return 1.0-nb0
+ return 1.0 - nb0
def calcul_ndvi(self):
"""
- Computer NDVI index for a Landsat image.
-
- NDVI = (PIR - R) / (PIR + R)
+ Function to compute NDVI index for a Landsat image with OTB BandMathX
+ - OTB help :
+ * il : Images list : data and cloud
+ * exp : Expression for ndvi : (PIR-R)/(PIR+R)
+ * out : Feature Output Image
"""
- np.seterr(invalid='ignore')
-
- # Extract raster's information
- data, in_ds = self.raster_data(self.cloudiness_pourcentage[0])
- data_cloud, in_ds_cloud = self.raster_data(self.cloudiness_pourcentage[1])
-
canal_PIR = Satellites.SATELLITE[self._class_archive._captor]["PIR"]
canal_R = Satellites.SATELLITE[self._class_archive._captor]["R"]
- # Computer NDVI
- mask = np.greater(data_cloud, -10000)
-
- # If True, -10000 (NaN) else compute mathematical operation
- ndvi = np.choose(mask, (-10000, eval('(data[canal_PIR]-data[canal_R])/(data[canal_PIR]+data[canal_R])')))
-
- # Outfile name
img_ndvi = "{0}_ndvi.TIF".format(self.cloudiness_pourcentage[0][:-4])
self.cloudiness_pourcentage.append(img_ndvi) # Add ndvi image path
- self.cloudiness_pourcentage.append(ndvi) # Add ndvi pixel value table
- self.create_empty_raster(img_ndvi, ndvi, in_ds)
+ otb_MathBand = otb.Registry.CreateApplication("BandMathX")
+ otb_MathBand.SetParameterStringList("il", [self.cloudiness_pourcentage[0], self.cloudiness_pourcentage[1]])
+ # otb_MathBand.SetParameterString("exp",\
+ # "(im1b1 > -1000 ? (im2b1 == 0 ? (im1b{0}-im1b{1})/(im1b{0}+im1b{1}) : -1000) : im1b1)"\
+ # .format(canal_PIR, canal_R))
+ otb_MathBand.SetParameterString("exp",\
+ "(im1b1 > -1 ? (im2b1 == 0 ? ndvi(im1b{1}, im1b{0}) : -1) : im1b1)"\
+ .format(canal_PIR, canal_R))
+ otb_MathBand.SetParameterString("out", img_ndvi)
+
+ otb_MathBand.ExecuteAndWriteOutput()
+
def create_empty_raster(self, out_raster_path, data, in_ds):
"""
@@ -330,8 +330,10 @@ class RasterSat_by_date(object):
int -- variable **nbband**, Band number of the out layer.
int -- variable **e**, Index to know if the raster exists.
- If it doesn't exists e = 0 else e = 1 (by default).
+ If it doesn't exists e = 0 else e = 1 (by default).
"""
+
+ # raise(BaseException("TEST"))
e = 1 # Raster out exists by default
@@ -349,6 +351,7 @@ class RasterSat_by_date(object):
if not os.path.exists(out_raster_path):
e = 0
+
# Create outfile
self.out_ds = driver.Create(out_raster_path, in_ds.RasterXSize, in_ds.RasterYSize, nbband, osgeo.gdal.GDT_Float32)
@@ -418,7 +421,7 @@ class RasterSat_by_date(object):
# Closing and statistics on output raster
out_band.FlushCache()
- out_band.SetNoDataValue(-10000)
+ # out_band.SetNoDataValue(-10000)
out_band.GetStatistics(-1, 1)
out_band = None
diff --git a/Sample.py b/Sample.py
index cb2d2f12d259bd6f029d8c00ce97afa3388d5f1c..b0ac80a79d00c382eee8a752b310c1b5859f69d8 100644
--- a/Sample.py
+++ b/Sample.py
@@ -82,6 +82,13 @@ class Sample(Vector):
if not os.path.exists(self.output_folder) :
os.mkdir(self.output_folder)
+
+ """
+ TODO : Enlever random.seed(0) pour retourner en aléatoire
+ """
+
+ random.seed(0)
+
# If the users want select polygons with a certain class name
if kw_field and kw_classes:
# The random sample with class name selected only
@@ -205,7 +212,7 @@ class Sample(Vector):
sys.exit(1)
# Specific output layer
- out_layer = out_ds.CreateLayer(str(output_sample), srsObj, geom_type=ogr.wkbMultiPolygon)
+ out_layer = out_ds.CreateLayer(output_sample, srsObj, geom_type=ogr.wkbMultiPolygon)
# Add existing fields
for i in range(0, len(self.field_names)):
@@ -229,7 +236,8 @@ class Sample(Vector):
in_feature = shp_ogr.SetNextByIndex(cnt)
in_feature = shp_ogr.GetNextFeature()
- geom = in_feature.GetGeometryRef() # Extract input geometry
+ # Extract input geometry
+ geom = in_feature.GetGeometryRef()
# Create a new polygon
out_feature = ogr.Feature(featureDefn)
@@ -237,11 +245,15 @@ class Sample(Vector):
# Set the polygon geometry and attribute
out_feature.SetGeometry(geom)
for i in range(0, len(self.field_names)):
- out_feature.SetField(self.field_names[i], in_feature.GetField(self.field_names[i]))
+ out_feature.SetField(self.field_names[i],\
+ float("%.5f" % in_feature.GetField(self.field_names[i])) \
+ if type(in_feature.GetField(self.field_names[i])) is float \
+ else in_feature.GetField(self.field_names[i]))
+
# In Option : Add a integer field
if add_field == 1:
out_feature.SetField(opt_field, opt_class[0])
-
+
# Append polygon to the output shapefile
out_layer.CreateFeature(out_feature)
diff --git a/Segmentation.py b/Segmentation.py
index de9a55507e205f55989117c7e4af9d5050341450..2ac1c1194fc68f616f48485b49633bbe743f6c7e 100644
--- a/Segmentation.py
+++ b/Segmentation.py
@@ -152,6 +152,7 @@ class Segmentation(Vector):
in_feature = shp_ogr.SetNextByIndex(0) # Polygons initialisation
in_feature = shp_ogr.GetNextFeature()
+ self.logger.debug(self.stats_rpg_tif)
#Â Loop on input polygons
while in_feature:
@@ -170,6 +171,7 @@ class Segmentation(Vector):
#Â Set the RPG column
recouv_crops_RPG = 0
pourc_inter = self.stats_rpg_tif[in_feature.GetFID()]['Maj_count_perc']
+
if pourc_inter >= 85:
recouv_crops_RPG = self.stats_rpg_tif[in_feature.GetFID()]['Maj_count']
@@ -292,22 +294,14 @@ class Segmentation(Vector):
distri_bio = []
distri_den = []
-
-
# for b in distri:
for idx, b in enumerate(distri):
if self.out_threshold["PREDICTION"][idx] in [1,2,8,9,10]\
- and b[-1] != -10000 and b[-1] < 1:
+ and b[-1] is not None and b[-1] != -10000 and b[-1] < 1:
distri_bio.append(b)
else:
distri_den.append(b)
-
- self.logger.debug(len(distri))
- self.logger.debug(len(distri_bio))
- self.logger.debug(len(distri_den))
- raise(BaseException("Error"))
-
#Â Transpose table
t_distri_bio = list(map(list, zip(*distri_bio)))
t_distri_den = list(map(list, zip(*distri_den)))
@@ -316,16 +310,21 @@ class Segmentation(Vector):
stdmore = (np.mean(t_distri_bio[2]) + np.std(t_distri_bio[2])) / np.max(t_distri_bio[2])
stdless = (np.mean(t_distri_bio[2]) - np.std(t_distri_bio[2])) / np.max(t_distri_bio[2])
- # self.out_threshold.append('<'+str(stdmore))
- # self.out_threshold.append('')
- # self.out_threshold.append('>'+str(stdless))
self.out_threshold["STDMORE"] = stdmore
self.out_threshold["STDLESS"] = stdless
# Compute density and biomass maximum
- self.max_wood_idm = np.max(t_distri_den[1])
- self.max_wood_sfs = np.max(t_distri_den[0])
- self.max_bio = np.max(t_distri_bio[2])
+ wood_sfs = np.array(t_distri_den[0], dtype=np.float64)
+ wood_idm = np.array(t_distri_den[1], dtype=np.float64)
+ max_bio = np.array(t_distri_den[2], dtype=np.float64)
+
+ self.max_wood_sfs = np.nanmax(wood_sfs)
+ self.max_wood_idm = np.nanmax(wood_idm)
+ self.max_bio = np.nanmax(max_bio)
+
+ # self.max_wood_sfs = np.max(t_distri_den[0])
+ # self.max_wood_idm = np.max(t_distri_den[1])
+ # self.max_bio = np.max(t_distri_bio[2])
def append_scale(self, select_class, form):
"""
@@ -342,7 +341,7 @@ class Segmentation(Vector):
for ind_stats in range(len(self.stats_dict)):
# Only valid on the second level
try:
- if self.class_tab_final[ind_stats][1] == select_class:
+ if self.class_tab_final[ind_stats][1] == select_class and self.stats_dict[ind_stats][3] is not None:
self.class_tab_final[ind_stats].insert(len(self.class_tab_final[ind_stats])-2,eval(form))
# To add phytomasse
try :
diff --git a/Toolbox.py b/Toolbox.py
index e5ccfc9012a239465e7f04dc6e4ef52534a2dbaf..8093cee1fb748c63bfe6043195bad2c14361e080 100644
--- a/Toolbox.py
+++ b/Toolbox.py
@@ -19,7 +19,9 @@
from datetime import date
-import os, sys, subprocess, glob
+import otbApplication as otb
+
+import os, sys, subprocess, glob, gc, psutil
import numpy as np
import json
import osgeo.gdal
@@ -77,10 +79,10 @@ class Toolbox(object):
return outclip
- def calc_serie_stats(self, table):
+ def calc_serie_stats(self, table, stats_name, output_folder):
"""
Function to compute stats on temporal cloud and ndvi spectral table
- Ndvi stats : min max std max-min
+ Ndvi stats : min max
@param table: Spectral data, cloud raster and ndvi raster
@type table: numpy.ndarray
@@ -90,68 +92,37 @@ class Toolbox(object):
numpy.ndarray -- variable **account_cloud**, pixel number clear on the area.
"""
- # Compute stats on these indexes
- ind = ['np.min(tab_ndvi_masked, axis=2)', 'np.max(tab_ndvi_masked, axis=2)', 'np.std(tab_ndvi_masked, axis=2)', \
- 'np.max(tab_ndvi_masked, axis=2)-np.min(tab_ndvi_masked, axis=2)'] # [Min, Max, Std, Max-Min]
+ ind = ['min({0})', 'max({0})']
+ no_data_value = [99999, -10000]
+
+ stack_ndvi = list(table[5])
+ image = "im{0}b1"
- # table = ['Année', 'Mois', 'Jour', chemin image', 'chemin image nuage',
- # 'données nuage', 'chemin image ndvi', 'données ndvi' ]
+ otb_MathBand = otb.Registry.CreateApplication("BandMathX")
+ otb_MathBand.SetParameterStringList("il", stack_ndvi)
- # For the cloud map
-
- # In the input table the cloud floor is the 6th
- tab_cloud = np.dstack(table[5]) # Stack cloud table (dimension : 12*X*Y to X*Y*12)
-
- # if tab_cloud = 0 then True else False / Mask cloud
- cloud_true = (tab_cloud == 0)
+ path_stats = []
- self.logger.debug('Mémoire cloud_true : {0}'.format(sys.getsizeof(cloud_true)))
+ for idx, i in enumerate(ind):
- # Account to tab_cloud if != 0 => Sum of True. (Dimension X*Y)
- account_cloud = np.sum(cloud_true, axis=2)
+ expression = "({0} != {1} ? {0} : im1b1)".format(\
+ i.format(",".join("({0}>-1000 ? {0} : {1})".format(\
+ image.format(j+1), no_data_value[idx]) for j in range(len(stack_ndvi)))), no_data_value[idx])
- cloud_true = None
+ out_ndvistats_raster = "{0}/{1}.TIF".format(output_folder, stats_name[idx])
- self.logger.debug('Mémoire cloud_true : {0}'.format(sys.getsizeof(cloud_true)))
-
- # For the ndvi stats
+ path_stats.append(out_ndvistats_raster)
- # In the input table the ndvi floor is the 8th
- stack_ndvi = np.dstack(table[7]) # Like cloud table, stack ndvi table
+ otb_MathBand.SetParameterString("exp",expression)
+ otb_MathBand.SetParameterString("out", out_ndvistats_raster)
- # Mask values -10000 and > 0 (Not cloud)
- mask_cloud = np.ma.masked_where((stack_ndvi == -10000) | (tab_cloud != 0), tab_cloud)
+ otb_MathBand.SetParameterOutputImagePixelType("out", otb.ImagePixelType_float)
- # Ndvi table with clear values
- tab_ndvi_masked = np.ma.array(stack_ndvi, mask=mask_cloud.mask)
-
- # Stats on the indexes defined above
- account_stats = []
+ otb_MathBand.ExecuteAndWriteOutput()
- for i in ind:
- i_stats = eval(i) # Compute stats
- i_stats.fill_value = -10000 # Substitute default fill value by -10000
- account_stats.append(i_stats.filled()) # Add stats table with true fill value
-
- #Â To extract index of the minimum ndvi to find the date
- if ind.index(i) == 0:
- i_date = eval(i.replace('np.min','np.argmin'))
-
- #Â To create a matrix with the date of pxl with a ndvi min
- tab_date = np.transpose(np.array(table[:3], dtype=object))
-
- # Loop on the temporal sequency
- for d in range(len(tab_date)):
- mask_data = np.ma.masked_values(i_date, d)
- i_date = mask_data.filled(date(int(tab_date[d][0]), int(tab_date[d][1]), int(tab_date[d][2])).toordinal()) # Date = day for year =1 day = 1 and month = 1
- mask_data = None
-
- account_stats.append(i_date) # Add the date table
+ return path_stats
- mask_cloud = None
- tab_ndvi_masked = None
- return account_stats, account_cloud
def check_proj(self):
"""
diff --git a/Vector.py b/Vector.py
index 4043d2151cfd60e10bb83a3aff6e2f2d9a072319..8725f7c125c5dee6f4933c008165c743c4e814db 100644
--- a/Vector.py
+++ b/Vector.py
@@ -21,6 +21,8 @@ import os, sys
import subprocess
import numpy as np
import osgeo
+from rpy2 import robjects
+import time
# try :
# import ogr, gdal
@@ -28,6 +30,7 @@ import osgeo
# from osgeo import ogr, gdal
import rasterstats
+import otbApplication as otb
from collections import *
import Outils
from RasterSat_by_date import RasterSat_by_date
@@ -138,17 +141,33 @@ class Vector():
:param (inraster,band): inraster -> Input image path, and band -> band number
:type (inraster,band): tuple
+
:kwargs: **rank** (int) - Zonal stats ranking launch
**nb_img** (int) - Number images launched with zonal stats
"""
+
ranking = kwargs['rank'] if kwargs.get('rank') else 0
nb_img = kwargs['nb_img'] if kwargs.get('nb_img') else 1
+ time1 = time.time()
self.logger.info('Compute {0} stats on {1} band {2}'.format(os.path.split(self.vector_used)[1], os.path.split(inraster)[1], band))
- stats = rasterstats.zonal_stats(self.vector_used, inraster, stats =['mean'], nodata=float('nan'), band=band)
+ zonal_stats = otb.Registry.CreateApplication("ZonalStatistics")
+ zonal_stats.SetParameterString("in", inraster)
+ zonal_stats.SetParameterString("vec", self.vector_used)
+ zonal_stats.SetParameterInt("band", band-1)
+ zonal_stats.SetParameterInt("ram", 128)
+ statistique = ["mean"]
+ zonal_stats.SetParameterStringList("stats", statistique)
+
+ zonal_stats.Execute()
+
+ stats = list(zonal_stats.GetParameterValue("out"))
+
+ # stats = rasterstats.zonal_stats(self.vector_used, inraster, stats =['mean'], nodata=float('nan'), band=band)
+
for i in range(len(stats)):
temp = defaultdict(lambda : [0]*nb_img)
for j in range(nb_img):
@@ -157,10 +176,15 @@ class Vector():
except IndexError:
pass
- temp[0][ranking] = list(stats)[i]['mean']
+ # temp[0][ranking] = list(stats)[i]['mean']
+ temp[0][ranking] = float(stats[i])
self.stats_dict[i] = temp[0]
-
- self.logger.info('End of stats on {0}'.format(os.path.split(inraster)[1]))
+ # self.logger.debug("self.stats_dict[i] : {0}".format(self.stats_dict[i]))
+
+ time2 = time.time()
+
+ self.logger.info('End of stats on {0} band {1}'.format(os.path.split(inraster)[1], band))
+ self.logger.debug("Temps d'éxecution {:.3f} ms".format((time2-time1)*1000.0))
def zonal_stats_pp(self, inraster):
"""
@@ -168,34 +192,76 @@ class Vector():
:param inraster: Input image path
:type inraster: str
+
:returns: dict -- **p_stats** : dictionnary with pxl percent in every polygon. Mainly 'Maj_count' (Majority Value) and 'Maj_count_perc' (Majority Percent)
"""
- p_stats = rasterstats.zonal_stats(self.vector_used, inraster, stats=['count'], copy_properties=True, categorical=True)
-
- for i in range(len(p_stats)):
- percent = 0.0
- for p in p_stats[i].keys():
- if type(p) == np.float32 and p_stats[i][p]/float(p_stats[i]['count'])*100 > percent:
- p_stats[i]['Maj_count'] = p
- p_stats[i]['Maj_count_perc'] = p_stats[i][p]/float(p_stats[i]['count'])*100
- percent = p_stats[i]['Maj_count_perc']
- if not 'Maj_count' in p_stats[i].keys() or not 'Maj_count_perc' in p_stats[i].keys():
- p_stats[i]['Maj_count']=0
- p_stats[i]['Maj_count_perc']=0
+ # p_stats = rasterstats.zonal_stats(self.vector_used, inraster, stats=['count'], copy_properties=True, categorical=True)
+
+ # for i in range(len(p_stats)):
+ # percent = 0.0
+ # for p in p_stats[i].keys():
+ # if type(p) == np.float32 and p_stats[i][p]/float(p_stats[i]['count'])*100 > percent:
+ # p_stats[i]['Maj_count'] = p
+ # p_stats[i]['Maj_count_perc'] = p_stats[i][p]/float(p_stats[i]['count'])*100
+ # percent = p_stats[i]['Maj_count_perc']
+
+ # if not 'Maj_count' in p_stats[i].keys() or not 'Maj_count_perc' in p_stats[i].keys():
+ # p_stats[i]['Maj_count']=0
+ # p_stats[i]['Maj_count_perc']=0
+
+ # for i in range(len(p_stats)):
+ # percent = 0.0
+ # Maj_count = 0
+ # Maj_count_perc = 0
+ # for p in p_stats[i].keys():
+ # if not isinstance(p, str) and p_stats[i][p]/float(p_stats[i]['count'])*100 > percent:
+ # Maj_count = p
+ # Maj_count_perc = p_stats[i][p]/float(p_stats[i]['count'])*100
+ # percent = Maj_count_perc
+
+ # p_stats[i]['Maj_count']= Maj_count
+ # p_stats[i]['Maj_count_perc']= Maj_count_perc
+
+ zonal_stats = otb.Registry.CreateApplication("ZonalStatistics")
+ zonal_stats.SetParameterString("in", inraster)
+ zonal_stats.SetParameterString("vec", self.vector_used)
+ zonal_stats.SetParameterInt("band", 1)
+ zonal_stats.SetParameterInt("ram", 128)
+ statistique = ["count"]
+ zonal_stats.SetParameterStringList("stats", statistique)
+
+ zonal_stats.Execute()
+
+ stats_otb = list(zonal_stats.GetParameterValue("out"))
+
+ p_stats = []
+
+ for item in stats_otb:
+ liste_item = item.split()
+ dico = dict()
+
+ dico['count'] = float(liste_item[0])
+ dico['Maj_count'] = float(liste_item[1])
+ dico['Maj_count_perc'] = float(liste_item[2])
+ p_stats.append(dico)
return p_stats
def layer_rasterization(self, raster_head, attribute_r, **kwargs):
"""
- Function to rasterize a vector. Define resolution, projection of the output raster with a raster head.
- And complete the gdal pointer empty properties with the layer's information of the vector and a defined field.
+ Function to rasterize a vector.
+ Define resolution, projection of the output raster with a raster head.
+ And complete the gdal pointer empty properties with the layer's information
+ of the vector and a defined field.
If a raster has several band, in option you can choice if you want one band or more.
:param raster_head: Raster path that will look like the final raster of the rasterization
:type raster_head: str
+
:param attribute_r: Field name of the shapefile that contains class names
:type attribute_r: str
+
:kwargs: **choice_nb_b** (int) - Output image number of band. If you choice 1, take first band. If you choice 2, take two first band etc...
:returns: str -- **valid_raster** : output raster path from the rasterization
"""
@@ -204,10 +270,12 @@ class Vector():
# for the validation shapefile
example_raster = RasterSat_by_date('', '') # Call the raster class
example_raster.choice_nb_b = kwargs['choice_nb_b'] if kwargs.get('choice_nb_b') else 0
+
raster_info = example_raster.raster_data(raster_head)# Extract data info
# Define the validation's vector
- valid_raster = self.vector_used[:-3]+'TIF' # Name of the output raster
+ valid_raster = "{0}.TIF".format(self.vector_used[:-4])# Name of the output raster
+
if os.path.exists(str(valid_raster)):
os.remove(valid_raster)
@@ -219,11 +287,12 @@ class Vector():
data_raster = raster_info[0]
info_out = example_raster.create_empty_raster(valid_raster, data_raster, raster_info[1])
+
self.raster_ds = example_raster.out_ds
# Virtual rasterize the vector
pt_rast = osgeo.gdal.RasterizeLayer(self.raster_ds, [1], self.data_source.GetLayer(), \
- options=["ATTRIBUTE=" + str(attribute_r)])
+ options=["ATTRIBUTE = {0}".format(attribute_r)])
if pt_rast != 0:
raise Exception("error rasterizing layer: %s" % pt_rast)
diff --git a/Vhrs.py b/Vhrs.py
index 99c3102dd8e2d9558c3be15a3442ed7c142fad88..7edd0f452adcc3d66feff45dd808cb5c8eb7ea3b 100644
--- a/Vhrs.py
+++ b/Vhrs.py
@@ -17,18 +17,17 @@
# You should have received a copy of the GNU General Public License
# along with PHYMOBAT 1.2. If not, see <http://www.gnu.org/licenses/>.
-import os
-import subprocess
+import os, sys
+import glob
import Outils
+import psutil
import otbApplication as otb
from multiprocessing import Process
-
class Vhrs():
-
"""
- Class to compute Haralick and SFS textures because of OTB application in command line
+ Class to compute Haralick and SFS textures with OTB
@param imag: The input image path to compute texture image
@type imag: str
@@ -47,36 +46,32 @@ class Vhrs():
Create a new 'Texture' instance
"""
- self._imag = imag
+ self._image = imag
self.mp = mp
self.logger = Outils.Log('Log', "VHRS")
+
+ self.out_sfs = "{0}_sfs.TIF".format(self._image[:-4])
+ self.out_haralick = "{0}_haralick.TIF".format(self._image[:-4])
+
+ ram_dispo = round((psutil.virtual_memory().available/pow(2,20))/2)
+
+ self.logger.debug("ram_dispo : {0}".format(ram_dispo))
- self.logger.info('SFS image')
- self.out_sfs = "{0}_sfs.TIF".format(self._imag[:-4])
-
- if not os.path.exists(self.out_sfs):
- self.logger.info("SFS image doesn't exist !")
- p_sfs = Process(target=self.sfs_texture_extraction)
- p_sfs.start()
+ p_sfs = Process(target=self.sfs_texture_extraction)
+ p_sfs.start()
- if mp == 0:
- p_sfs.join()
-
- self.logger.info('Haralick image')
- self.out_haralick = "{0}_haralick.TIF".format(self._imag[:-4])
+ if mp == 0:
+ p_sfs.join()
- if not os.path.exists(self.out_haralick):
- self.logger.info("Haralick image doesn't exist !")
- p_har = Process(target=self.haralick_texture_extraction, args=('simple', ))
- p_har.start()
+ p_har = Process(target=self.haralick_texture_extraction, args=('simple', ))
+ p_har.start()
- if mp == 0:
- p_har.join()
+ if mp == 0:
+ p_har.join()
if mp == 1:
- if not os.path.exists(self.out_sfs) and not os.path.exists(self.out_haralick):
- p_sfs.join()
- p_har.join()
+ p_sfs.join()
+ p_har.join()
def sfs_texture_extraction(self):
"""
@@ -92,18 +87,24 @@ class Vhrs():
- parameters.maxcons : Ratio Maximum Consideration Number
* out : Feature Output Image
"""
+
+ self.logger.info('SFS image')
+
+ if not os.path.exists(self.out_sfs):
+
+ self.logger.info("SFS image doesn't exist !")
- otb_SFS = otb.Registry.CreateApplication("SFSTextureExtraction")
- otb_SFS.SetParameterString("in", self._imag)
- otb_SFS.SetParameterInt("channel", 2)
- otb_SFS.SetParameterInt("parameters.spethre", 50)
- otb_SFS.SetParameterInt("parameters.spathre", 100)
- otb_SFS.SetParameterString("out", self.out_sfs)
+ otb_SFS = otb.Registry.CreateApplication("SFSTextureExtraction")
+ otb_SFS.SetParameterString("in", self._image)
+ otb_SFS.SetParameterInt("channel", 2)
+ otb_SFS.SetParameterInt("parameters.spethre", 50)
+ otb_SFS.SetParameterInt("parameters.spathre", 100)
+ otb_SFS.SetParameterString("out", self.out_sfs)
- otb_SFS.ExecuteAndWriteOutput()
+ otb_SFS.ExecuteAndWriteOutput()
+
+ self.logger.info('SFS image created.')
- self.logger.info('SFS image created.')
-
def haralick_texture_extraction(self, texture_choice):
"""
Function to compute OTB Haralick texture image
@@ -129,15 +130,20 @@ class Vhrs():
@param texture_choice: Order texture choice -> Simple / Advanced / Higher
@type texture_choice: str
"""
-
- otb_Haralick = otb.Registry.CreateApplication("HaralickTextureExtraction")
- otb_Haralick.SetParameterString("in", self._imag)
- otb_Haralick.SetParameterInt("channel", 2)
- otb_Haralick.SetParameterInt("parameters.xrad", 3)
- otb_Haralick.SetParameterInt("parameters.yrad", 3)
- otb_Haralick.SetParameterString("texture", texture_choice)
- otb_Haralick.SetParameterString("out", self.out_haralick)
+ self.logger.info('Haralick image')
+
+ if not os.path.exists(self.out_haralick):
+
+ self.logger.info("Haralick image doesn't exist !")
+
+ otb_Haralick = otb.Registry.CreateApplication("HaralickTextureExtraction")
+ otb_Haralick.SetParameterString("in", self._image)
+ otb_Haralick.SetParameterInt("channel", 2)
+ otb_Haralick.SetParameterInt("parameters.xrad", 3)
+ otb_Haralick.SetParameterInt("parameters.yrad", 3)
+ otb_Haralick.SetParameterString("texture", texture_choice)
+ otb_Haralick.SetParameterString("out", self.out_haralick)
- otb_Haralick.ExecuteAndWriteOutput()
+ otb_Haralick.ExecuteAndWriteOutput()
- self.logger.info('Haralick image created.')
\ No newline at end of file
+ self.logger.info('Haralick image created.')