Code source de Processing

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# This file is part of PHYMOBAT 2.0.
# Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS)
# 
# PHYMOBAT 2.0 is free software: you can redistribute it and/or modify
# 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 2.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 2.0.  If not, see <http://www.gnu.org/licenses/>.

import os, sys, math
import numpy as np
import subprocess
from sklearn.ensemble import RandomForestClassifier
from sklearn import tree
try :
    import ogr, gdal
except :
    from osgeo import ogr, gdal

from Toolbox import Toolbox
from Seath import Seath
from Precision_moba import Precision_moba
# Class group image
from Archive import Archive
from RasterSat_by_date import RasterSat_by_date
from Vhrs import Vhrs
from Slope import Slope

# Class group vector
from Vector import Vector
from Sample import Sample
from Segmentation import Segmentation
from Rpg import Rpg

from collections import defaultdict
from multiprocessing import Process
from multiprocessing.managers import BaseManager, DictProxy

[docs]class Processing(): """ Main processing. This class launch the others system classes. It take into account CarHab classification method MOBA. This way is broken down into 3 parts : - Image Processing (Search, download and processing) - Vector Processing (Optimal threshold, Sample processing) - Classification - Validation **Main parameters** :param captor_project: Satellite captor name :type captor_project: str :param classif_year: Classification year :type classif_year: str :param nb_avalaible_images: Number download available images :type nb_avalaible_images: int :param path_folder_dpt: Main folder path :type path_folder_dpt: str :param folder_archive: Archive downloaded folder path :type folder_archive: str :param folder_processing: Processing folder name. By default : 'Traitement' :type folder_processing: str :param path_area: Study area shapefile :type path_area: str :param path_ortho: VHRS image path :type path_ortho: str :param path_mnt: MNT image path :type path_mnt: str :param path_segm: Segmentation shapefile :type path_segm: str **Id information to download on theia platform** :param user: Connexion Username :type user: str :param password: Connexion Password :type password: str **Output parameters** :param output_name_moba: Output classification shapefile :type output_name_moba: str :param out_fieldname_carto: Output shapefile field name :type out_fieldname_carto: list of str :param out_fieldtype_carto: Output shapefile field type :type out_fieldtype_carto: list of str (eval ogr pointer) **Sample parameters** :param fieldname_args: Sample field names 2 by 2 :type fieldname_args: list of str :param class_args: Sample class names 2 by 2 :type class_args: list of str :param sample_name: List of sample name (path) :type sample_name: list of str :param list_nb_sample: Number of polygons for every sample :type list_nb_sample: list of int **Multi-processing parameters** :param mp: Boolean variable -> 0 or 1. - 0 means, not multi-processing - 1 means, launch process with multi-processing :type mp: int """ def __init__(self): # Used variables self.captor_project = '' self.classif_year = '' self.path_folder_dpt = '' self.folder_archive = '' self.folder_processing = 'Traitement' self.path_area = '' self.path_ortho = '' self.path_mnt = '' self.path_segm = '' self.output_name_moba = '' self.w_proxy = None # For the "Proxy window" # Id information to download on theia platform self.user = '' self.password = '' # List of output raster path self.raster_path = [] self.list_band_outraster = [] # Class name # TODO : Change index of the classes -> Harbacées 6 / Ligneux 7 by Agriculuture 4 / Eboulis 5 self.in_class_name = ['Non Vegetation semi-naturelle', 'Vegetation semi-naturelle',\ 'Herbacees', 'Ligneux', \ 'Ligneux mixtes', 'Ligneux denses',\ 'Agriculture', 'Eboulis', \ 'Forte phytomasse', 'Moyenne phytomasse', 'Faible phytomasse'] # Sample field names 2 by 2 self.fieldname_args = [] # 'CODE_GROUP', 'CODE_GROUP',\ # 'echant', 'echant',\ # 'echant', 'echant'] # Sample class names 2 by 2 self.class_args = [] # '1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 14, 15, 20, 21, 22, 23, 24, 26, 28', '18',\ # 'H','LF, LO',\ # 'LF', 'LO'] # Decision tree combination self.tree_direction = [[0],\ [0],\ [1, 3, 4],\ [1, 3, 5],\ [1, 2, 8],\ [1, 2, 9],\ [1, 2, 10]] # [['Cultures'],['Semi-naturelles', 'Herbacees', 'Forte phytomasse'], ... # ..., ['Semi-naturelles', 'Ligneux', 'Ligneux denses']] # Slope degrees self.slope_degree = 30 # Output shapefile field name self.out_fieldname_carto = ['ID', 'AREA'] #, 'NIVEAU_1', 'NIVEAU_2', 'NIVEAU_3', 'POURC'] # Output shapefile field type self.out_fieldtype_carto = [ogr.OFTString, ogr.OFTReal] #, ogr.OFTString, ogr.OFTString, ogr.OFTString, ogr.OFTReal] # List of sample name path self.sample_name = [] # List of number sample self.list_nb_sample = [] # Number download available images self.nb_avalaible_images = 0 # Multi-processing variable self.mp = 1 # Function followed self.check_download = '' self.decis = {} # Images after processing images self.out_ndvistats_folder_tab = defaultdict(list) # Validation shapefiles information self.valid_shp = [] # Radom Forest Model # Set the parameters of this random forest from the estimator self.rf = RandomForestClassifier(n_estimators=500, criterion='gini', max_depth=None, min_samples_split=2, \ min_samples_leaf=1, max_features='auto', \ bootstrap=True, oob_score=True)
[docs] def i_tree_direction(self): """ Interface function to can extract one level or two levels of the final classification """ if len(self.out_fieldname_carto) == 3: self.tree_direction = [[0], [1]] if len(self.out_fieldname_carto) == 4: self.tree_direction = [[0], [0], [1, 2], [1, 3]]
[docs] def i_download(self, dd): """ Interface function to download archives on the website Theia Land. This function extract the number of downloadable image with :func:`Archive.Archive.listing`. Then, this function download :func:`Archive.Archive.download` and unzip :func:`Archive.Archive.decompress` images in the archive folder (**folder_archive**). :param dd: Boolean variable to launch download images -> 0 or 1. - 0 means, not downloading - 1 means, launch downloading :type dd: int """ self.folder_archive = self.captor_project + '_PoleTheia' self.check_download = Archive(self.captor_project, self.classif_year, self.path_area, self.path_folder_dpt, self.folder_archive, self.w_proxy) self.check_download.listing() self.nb_avalaible_images = len(self.check_download.list_archive) # check_download.set_list_archive_to_try(check_download.list_archive[:3]) if dd == 1: # self.check_download.download_auto(self.user, self.password) self.check_download.download_auto(self.user, self.password) self.check_download.decompress()
[docs] def i_glob(self): """ Function to load existing images to launch the processing. It need to archives. Then, to select processing images, select archives """ self.folder_archive = self.captor_project + '_PoleTheia' self.check_download = Archive(self.captor_project, self.classif_year, self.path_area, self.path_folder_dpt, self.folder_archive, self.w_proxy) self.check_download.decompress()
[docs] def i_img_sat(self): """ Interface function to processing satellite images: 1. Clip archive images and modify Archive class to integrate clip image path. With :func:`Toolbox.clip_raster` in ``Toolbox`` module. 2. Search cloud's percentage :func:`RasterSat_by_date.RasterSat_by_date.pourc_cloud`, select image and compute ndvi index :func:`RasterSat_by_date.RasterSat_by_date.calcul_ndvi`. If cloud's percentage is greater than 40%, then not select and compute ndvi index. 3. Compute temporal stats on ndvi index [min, max, std, min-max]. With :func:`Toolbox.calc_serie_stats` in ``Toolbox`` module. 4. Create stats ndvi raster and stats cloud raster. >>> import RasterSat_by_date >>> stats_test = RasterSat_by_date(class_archive, big_folder, one_date) >>> stats_test.complete_raster(stats_test.create_raster(in_raster, stats_data, in_ds), stats_data) """ # Clip archive images and modify Archive class to integrate clip image path for clip in self.check_download.list_img: clip_index = self.check_download.list_img.index(clip) current_list = Toolbox() current_list.imag = clip[3] current_list.vect = self.path_area current_list.check_proj() # Check if projection is RFG93 self.check_download.list_img[clip_index][3] = current_list.clip_raster() # Multispectral images current_list.imag = clip[4] current_list.check_proj() # Check if projection is RFG93 self.check_download.list_img[clip_index][4] = current_list.clip_raster() # Cloud images # Images pre-processing spectral_out = [] for date in self.check_download.single_date: check_L8 = RasterSat_by_date(self.check_download, self.folder_processing, date) check_L8.mosaic_by_date() # Search cloud's percentage, select image and compute ndvi index if > cl cl = check_L8.pourc_cloud(check_L8._one_date[3], check_L8._one_date[4]) if cl > 0.60: check_L8.calcul_ndvi(check_L8._one_date[3]) spectral_out.append(check_L8._one_date) # 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) # Create stats ndvi raster and stats cloud raster stats_L8 = RasterSat_by_date(self.check_download, self.folder_processing, [0]) # Stats cloud raster out_cloud_folder = stats_L8._class_archive._folder + '/' + stats_L8._big_folder + '/' + self.captor_project + \ '/Cloud_number_' + self.captor_project + '.TIF' stats_L8.complete_raster(stats_L8.create_raster(out_cloud_folder, stats_cloud, \ stats_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_folder = stats_L8._class_archive._folder + '/' + stats_L8._big_folder + '/' + self.captor_project + \ '/' + stats_name[stats_index] + '_' + self.captor_project + '.TIF' self.out_ndvistats_folder_tab[stats_index] = out_ndvistats_folder stats_L8.complete_raster(stats_L8.create_raster(out_ndvistats_folder, stats_ndvi[stats_index], \ stats_L8.raster_data(self.check_download.list_img[0][4])[1]), \ stats_ndvi[stats_index])
[docs] def i_slope(self): """ Interface function to processing slope raster. From a MNT, and with a command line gdal, this function compute slope in degrees :func:`Slope.Slope`. """ current_path_mnt = Toolbox() current_path_mnt.imag = self.path_mnt current_path_mnt.vect = self.path_area path_mnt = current_path_mnt.clip_raster() study_slope = Slope(path_mnt) study_slope.extract_slope()# Call this function to compute slope raster self.path_mnt = study_slope.out_mnt
[docs] def i_vhrs(self):#, vs): """ Interface function to processing VHRS images. It create two OTB texture images :func:`Vhrs.Vhrs` : SFS Texture and Haralick Texture """ # Create texture image # Clip orthography image current_path_ortho = Toolbox() current_path_ortho.imag = self.path_ortho current_path_ortho.vect = self.path_area 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['haralick'] = texture_irc.out_haralick
[docs] def i_images_processing(self, vs): """ Interface function to launch processing VHRS images :func:`i_vhrs` and satellite images :func:`i_img_sat` in multi-processing. :param vs: Boolean variable to launch processing because of interface checkbox -> 0 or 1. - 0 means, not texture processing - 1 means, launch texture processing :type vs: int """ # Multiprocessing mgr = BaseManager() mgr.register('defaultdict', defaultdict, DictProxy) mgr.start() self.out_ndvistats_folder_tab = mgr.defaultdict(list) p_img_sat = Process(target=self.i_img_sat) p_img_sat.start() if self.mp == 0: p_img_sat.join() if vs == 1: p_vhrs = Process(target=self.i_vhrs)#, args=(vs, )) p_vhrs.start() p_vhrs.join() if self.mp == 1: p_img_sat.join() # List of output raster path self.raster_path.append(self.out_ndvistats_folder_tab[0]) # List of output raster band self.list_band_outraster.append(1) if vs == 1: self.raster_path.append(self.out_ndvistats_folder_tab['sfs']) self.list_band_outraster.append(4) self.raster_path.append(self.out_ndvistats_folder_tab['haralick']) self.list_band_outraster.append(2) # To slope, to extract scree if self.path_mnt != '': self.raster_path.append(self.path_mnt) self.list_band_outraster.append(1) self.raster_path.append(self.out_ndvistats_folder_tab[2]) # example raster path tab : # [path_folder_dpt + '/' + folder_processing + '/' + classif_year + '/Min_2014.TIF',\ # os.path.dirname(path_ortho) + '/Clip_buffer_surface_dep_18_IRCOrtho65_2m_sfs.TIF',\ # os.path.dirname(path_ortho) + '/Clip_buffer_surface_dep_18_IRCOrtho65_2m_haralick.TIF',\ # path_folder_dpt + '/' + folder_processing + '/' + classif_year + '/Max_2014.TIF'] # List of output raster band self.list_band_outraster.append(1) #[1, 4, 2, 1] print("End of images processing !")
[docs] def i_rpg(self, path_rpg): """ Interface function to extract mono rpg crops. :param path_rpg: Input RPG shapefile. :type path_rpg: str :returns: str -- variable **Rpg.vector_used**, output no duplicated crops shapefile (path). """ # Extract mono rpg crops mono_sample = Rpg(path_rpg, self.path_area) # If exists, do not create a mono rpg file if os.path.basename(str(path_rpg))[:5]!='MONO_': mono_sample.mono_rpg() mono_sample.create_new_rpg_files() else: print('MONO RPG file exists already !!!') print('End of RPG processing') return mono_sample.vector_used
[docs] def i_sample(self): """ Interface function to compute threshold with various sample. It also extract a list of validation layer (shapefile) to compute the precision of the next classification :func:`i_validate`. It create samples 2 by 2 with kwargs field names and class :func:`Sample.Sample.create_sample`. Then, it compute zonal statistics by polygons :func:`Vector.Sample.zonal_stats`. With zonal statistics computed, a optimal threshold is determined :func:`Seath.Seath.separability_and_threshold` that will print in a text file .lg in the main folder. .. warning:: :func:`Seath.Seath.separability_and_threshold` does not always allow to discriminate optimal threshold. Then, this function will be launch at least ten time until it reaches a optimal threshold. """ # Compute threshold with various sample i_s = 0 while i_s < 10: try : self.valid_shp = [] sample_rd = {} for sple in range(len(self.sample_name) * 2): kwargs = {} kwargs['fieldname'] = self.fieldname_args[sple] kwargs['class'] = self.class_args[sple] sample_rd[sple] = Sample(self.sample_name[sple/2], self.path_area, self.list_nb_sample[sple/2]) sample_rd[sple].create_sample(**kwargs) sample_rd[sple].zonal_stats((self.raster_path[sple/2], self.list_band_outraster[sple/2])) # Add the validation shapefile self.valid_shp.append([sample_rd[sple].vector_val, kwargs['fieldname'], kwargs['class']]) # Search the optimal threshold by class # Open a text file to print stats of Seath method file_J = self.path_folder_dpt + '/log_J.lg' f = open(file_J, "wb") for th_seath in range(len(self.sample_name)): self.decis[th_seath] = Seath() self.decis[th_seath].value_1 = sample_rd[th_seath*2].stats_dict self.decis[th_seath].value_2 = sample_rd[th_seath*2 + 1].stats_dict self.decis[th_seath].separability_and_threshold() # Print the J value in the text file .lg f.write('For ' + str(self.sample_name[th_seath]) + ' :\n') f.write('J = ' + str(self.decis[th_seath].J[0]) +'\n') f.write('The class 1 ' + str(self.decis[th_seath].threshold[0]) +'\n') f.close() i_s = 20 except: i_s = i_s + 1 # Method to stop the processus if there is not found a valid threshold if i_s != 20: print 'Problem in the sample processing !!!' sys.exit(1)
[docs] def i_sample_rf(self): """ This function build a random forest trees like model to create a final classification. All of This using the method described in the :func:`i_validate` function and because of sklearn module. """ X_rf = [] y_rf = [] sample_rd = {} # Tricks to add all textural indexes rm_index = 1 self.raster_path.remove(self.raster_path[rm_index]) # Remove SFS layer self.raster_path.remove(self.raster_path[rm_index]) # Remove Haralick layer self.list_band_outraster.remove(self.list_band_outraster[rm_index]) # Remove band of the layer self.list_band_outraster.remove(self.list_band_outraster[rm_index]) # Remove band of the layer # Add all layers in the simple index haralick for add_layer in range(8): self.raster_path.insert(add_layer+1, self.out_ndvistats_folder_tab['haralick']) self.list_band_outraster.insert(add_layer+1, add_layer+1) # Add all layers in the SFS index for add_layer in range(6): self.raster_path.insert(add_layer+1, self.out_ndvistats_folder_tab['sfs']) self.list_band_outraster.insert(add_layer+1, add_layer+1) # Extract value mean from polygons for sple in range(len(self.sample_name) * 2): kwargs = {} kwargs['fieldname'] = self.fieldname_args[sple] kwargs['class'] = self.class_args[sple] sample_rd[sple] = Sample(self.sample_name[sple/2], self.path_area, self.list_nb_sample[sple/2]) sample_rd[sple].create_sample(**kwargs) # Add the validation shapefile self.valid_shp.append([sample_rd[sple].vector_val, kwargs['fieldname'], kwargs['class']]) for lbo in range(len(self.raster_path)): kwargs['rank'] = lbo kwargs['nb_img'] = len(self.raster_path) sample_rd[sple].zonal_stats((self.raster_path[lbo], self.list_band_outraster[lbo]), **kwargs) # To convert the dictionnary in a list for key, value in sample_rd[sple].stats_dict.iteritems(): # 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: # y_rf.append(1) pass elif sple == 3: # y_rf.append(4) pass else: y_rf.append(sple) X_rf.append([-10000 if (math.isnan(x) or math.isinf(x)) else x for x in value]) # Build a forest of trees from the samples self.rf = self.rf.fit(X_rf, y_rf) # Print in a file feature important importance = self.rf.feature_importances_ importance = [(importance[x],x+1) for x in range(len(importance))] importance.sort() file_feat_import = os.path.dirname(str(self.raster_path[0])) + '/Feature_important_RF.ft' if os.path.exists(file_feat_import): os.remove(file_feat_import) f_out = open(file_feat_import, "wb") f_out.write(str(importance)) # Close the output file f_out.close() # Print in a file decision tree file_decisiontree = os.path.dirname(str(self.raster_path[0])) + '/Decision_tree.dot' if os.path.exists(file_decisiontree): os.remove(file_decisiontree) tree_in_forest = self.rf.estimators_[499] with open(file_decisiontree, 'w') as my_file: my_file = tree.export_graphviz(tree_in_forest, out_file = my_file)
[docs] def i_classifier_rf(self): """ Interface function to launch random forest classification with a input segmentation :func:`Segmentation.Segmentation`. This function use the sklearn module to build the best of decision tree to extract classes. The optimal threshold are stored by class **rf** variable in :func:`Processing.i_sample_rf`. Then it computes zonal statistics by polygons for every images in multi-processing (if **mp** = 1). """ # Multiprocessing mgr = BaseManager() mgr.register('defaultdict', defaultdict, DictProxy) mgr.start() multi_process_var = [] # Multi processing variable # Extract final cartography out_carto = Segmentation(self.path_segm, self.path_area) out_carto.output_file = self.output_name_moba out_carto.out_class_name = self.in_class_name # out_carto.out_threshold = [] for ind_th in range(len(self.raster_path)): multi_process_var.append([self.raster_path[ind_th], self.list_band_outraster[ind_th]]) # Compute zonal stats with multi processing exist_stats = 1 # By default, the stats file exists already file_stats = os.path.dirname(str(self.raster_path[0])) + '/Stats_raster_spectral_texture.stats' # Stats backup file if not os.path.exists(file_stats): 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, "wb") p = [] kwargs = {} X_out_rf = [] # Variable list to compute decision tree with random forest method if exist_stats == 0: out_carto.stats_dict = mgr.defaultdict(list) for i in range(len(multi_process_var)): kwargs['rank'] = i kwargs['nb_img'] = len(multi_process_var) p.append(Process(target=out_carto.zonal_stats, args=(multi_process_var[i], ), kwargs=kwargs)) p[i].start() if self.mp == 0: p[i].join() if self.mp == 1: for i in range(len(multi_process_var)): p[i].join() for key, value_seg in out_carto.stats_dict.items(): true_value = [-10000 if (math.isnan(x) or math.isinf(x)) else x for x in value_seg] X_out_rf.append(true_value) # Print rasters stats value in the text file .lg f_out.write(str(true_value) + '\n') # Close the output file f_out.close() else: # If the stats file exists already, open this file and append in the stats_dict variable out_carto.stats_dict = defaultdict(list) with open(file_stats, "r") as f_in: index_in_stats=-1 for x_in in f_in.readlines(): 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) # For the higher than level 1 if len(self.sample_name) > 2: # Compute the biomass and density distribution # Use 'out_carto.out_threshold' to konw predicted in the segmentation class out_carto.out_threshold = predicted_rf # In the compute_biomass_density function, this variable used normally to define # threshold of the classification with SEATH method is initialized out_carto.compute_biomass_density('RF') out_carto.class_tab_final = defaultdict(list) for i_polyg in range(len(predicted_rf)): i_final = 0 class_final = [] # Initialize the predicted output format # For example : predicted => 4, formatted => [1,3,4] while i_final < len(self.tree_direction): if self.tree_direction[i_final][len(self.tree_direction[i_final])-1] == predicted_rf[i_polyg]: class_final = self.tree_direction[i_final] i_final = len(self.tree_direction) i_final = i_final + 1 if class_final == []: class_final = [1, 2] # Set the class name because of predicted output formatted out_carto.class_tab_final[i_polyg] = [self.in_class_name[f] for f in class_final] + \ [predicted_rf[i_polyg]] + [predicted_rf[i_polyg]] # For the output line with one level, add a phantom level # TODO : Replace constant by a variable in the condition 'while' while len(out_carto.class_tab_final[i_polyg]) < 5: out_carto.class_tab_final[i_polyg].insert(len(out_carto.class_tab_final[i_polyg])-2,'') # If there is more one fieldnames line edit fulled in classification tab if len(self.sample_name) > 2: # Compute biomass and density scale out_carto.append_scale(self.in_class_name[2], 'self.stats_dict[ind_stats][3]/self.max_bio') out_carto.append_scale(self.in_class_name[3], 'self.stats_dict[ind_stats][2]/self.max_wood_idm') # Rasterize RPG shapefile to complete the final shapefile opt = {} opt['Remove'] = 1 rpg_tif = Vector(self.sample_name[0], self.path_area, **opt) # if not os.path.exists(str(rpg_tif.vector_used[:-3]+'TIF')): kwargs['choice_nb_b'] = 1 out_carto.stats_rpg_tif = out_carto.zonal_stats_pp(rpg_tif.layer_rasterization(self.path_ortho, 'CODE_GROUP', **kwargs)) # Final cartography out_carto.create_cartography(self.out_fieldname_carto, self.out_fieldtype_carto)
[docs] def i_classifier_s(self): """ Interface function to launch decision tree classification with a input segmentation :func:`Segmentation.Segmentation`. This function store optimal threshold by class **Segmentation.out_threshold**. Then it computes zonal statistics by polygons for every images in multi-processing (if **mp** = 1). """ # Multiprocessing mgr = BaseManager() mgr.register('defaultdict', defaultdict, DictProxy) mgr.start() multi_process_var = [] # Multi processing variable # Extract final cartography out_carto = Segmentation(self.path_segm, self.path_area) out_carto.output_file = self.output_name_moba out_carto.out_class_name = self.in_class_name out_carto.out_threshold = [] for ind_th in range(len(self.sample_name)): out_carto.out_threshold.append(self.decis[ind_th].threshold[0]) if '>' in self.decis[ind_th].threshold[0]: out_carto.out_threshold.append(self.decis[ind_th].threshold[0].replace('>', '<=')) elif '<' in self.decis[ind_th].threshold[0]: out_carto.out_threshold.append(self.decis[ind_th].threshold[0].replace('<', '>=')) # out_carto.zonal_stats((raster_path[ind_th], list_band_outraster[ind_th])) multi_process_var.append([self.raster_path[ind_th], self.list_band_outraster[ind_th]]) # Compute zonal stats on slope raster multi_process_var.append([self.raster_path[ind_th+1], self.list_band_outraster[ind_th+1]]) out_carto.out_threshold.append('<'+str(self.slope_degree)) # To agriculture out_carto.out_threshold.append('>='+str(self.slope_degree)) # To scree if self.path_mnt != '': # Add class indexes self.tree_direction[0].append(6) self.tree_direction[0].append(7) # Compute zonal stats on Max NDVI raster try: # out_carto.zonal_stats((raster_path[ind_th+1], list_band_outraster[ind_th+1])) multi_process_var.append([self.raster_path[ind_th+2], self.list_band_outraster[ind_th+2]]) # Compute stats twice, because there is 3 classes and not 2 # out_carto.zonal_stats((raster_path[ind_th+1], list_band_outraster[ind_th+1])) multi_process_var.append([self.raster_path[ind_th+2], self.list_band_outraster[ind_th+2]]) except: print('Not MNT on the 3rd step') multi_process_var.append([self.raster_path[ind_th+1], self.list_band_outraster[ind_th+1]]) multi_process_var.append([self.raster_path[ind_th+1], self.list_band_outraster[ind_th+1]]) # Compute zonal stats with multi processing exist_stats = 1 # By default, the stats file exists already file_stats = os.path.dirname(str(self.raster_path[0])) + '/Stats_raster_spectral_texture.stats' # Stats backup file if not os.path.exists(file_stats): 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, "wb") p = [] kwargs = {} X_out_rf = [] # Variable list to compute decision tree with random forest method if exist_stats == 0: out_carto.stats_dict = mgr.defaultdict(list) for i in range(len(multi_process_var)): kwargs['rank'] = i kwargs['nb_img'] = len(multi_process_var) p.append(Process(target=out_carto.zonal_stats, args=(multi_process_var[i], ), kwargs=kwargs)) p[i].start() if self.mp == 0: p[i].join() if self.mp == 1: for i in range(len(multi_process_var)): p[i].join() for key, value_seg in out_carto.stats_dict.items(): true_value = [-10000 if (math.isnan(x) or math.isinf(x)) else x for x in value_seg] # Print rasters stats value in the text file .lg f_out.write(str(true_value) + '\n') # Close the output file f_out.close() else: # If the stats file exists already, open this file and append in the stats_dict variable out_carto.stats_dict = defaultdict(list) with open(file_stats, "r") as f_in: index_in_stats=-1 for x_in in f_in.readlines(): 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'))) # For the higher than level 1 if len(self.sample_name) > 2: # Compute the biomass and density distribution out_carto.compute_biomass_density() out_carto.class_tab_final = defaultdict(list) self.i_tree_direction() out_carto.decision_tree(self.tree_direction) # If there is more one fieldnames line edit fulled in classification tab if len(self.sample_name) > 2: # Compute biomass and density scale out_carto.append_scale(self.in_class_name[2], 'self.stats_dict[ind_stats][3]/self.max_bio') out_carto.append_scale(self.in_class_name[3], 'self.stats_dict[ind_stats][2]/self.max_wood_idm') # Rasterize RPG shapefile to complete the final shapefile opt = {} opt['Remove'] = 1 rpg_tif = Vector(self.sample_name[0], self.path_area, **opt) rpg_tif.layer_rasterization(self.path_ortho, 'CODE_GROUP') # Final cartography out_carto.mono_rpg_tif = self.sample_name[0][:-4] + '.TIF' out_carto.create_cartography(self.out_fieldname_carto, self.out_fieldtype_carto)
[docs] def i_validate(self): """ Interface to validate a classification. It going to rasterize the validation shapefile and the classification shapefile with :func:`layer_rasterization`. Next, to compare pixel by pixel, the classification quality to built a confusion matrix in a csv file. """ # Variable to convert the input classname to an individual interger # Only for the validate sample class_validate = 0 complete_validate_shp = os.path.dirname(str(self.valid_shp[0][0])) + '/validate.shp' # TODO: Set this method in the Precision_moba class # Processing to rasterize the validate shapefile. 1) Merge sahpefiles 2) Rasterization for val in self.valid_shp: if class_validate != 2: # Grassland to 1 if (class_validate !=3 and len(self.out_fieldname_carto) != 4+2) or len(self.out_fieldname_carto) == 4+2: # To the level 3 with woodeen to 4 and 5 # # Self.valid_shp is a list of list. In this variable there is : # [Shapefile path, fieldname classes, classnames] opt = {} opt['Remove'] = 1 # To overwrite # Create a raster to valide the classification # First time, create a new shapefile with a new field integer sample_val = Sample(val[0], self.path_area, 1, **opt) opt['add_fieldname'] = 1 opt['fieldname'] = 'CLASS_CODE' opt['class'] = str(class_validate) # Add integer classes # Set the new shapefile val[0] = val[0][:-4] + '_.shp' val[1] = opt['fieldname'] val[2] = opt['class'] # Complete the new shapefile sample_val.fill_sample(val[0], 0, **opt) # Second time, merge the validate shapefile if class_validate == 0: process_tocall_merge = ['ogr2ogr', '-overwrite', complete_validate_shp, val[0]] elif class_validate > 0: process_tocall_merge = ['ogr2ogr', '-update', '-append', complete_validate_shp, \ val[0], '-nln', os.path.basename(str(complete_validate_shp[:-4]))] subprocess.call(process_tocall_merge) # Increrment variable class_validate = self.valid_shp.index(val) + 1 # Compute precision of the classification valid = Precision_moba(self.path_area, self.path_folder_dpt) valid.complete_validation_shp = complete_validate_shp valid.ex_raster = self.raster_path[0] # TODO: Call the RasterSat_by_Date class here instead of the Precision_moba class valid.preprocess_to_raster_precision(self.output_name_moba, 'FBPHY_CODE') # To the classification's data valid.preprocess_to_raster_precision(complete_validate_shp, val[1]) # To the validation's data # Compute precision on the output classification valid.confus_matrix(valid.complete_img[0].raster_data(valid.img_pr[0])[0], \ valid.complete_img[1].raster_data(valid.img_pr[1])[0])