Code source de Segmentation

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# This file is part of PHYMOBAT 1.2.
# Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS)
# 
# PHYMOBAT 1.2 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 1.2 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 1.2.  If not, see <http://www.gnu.org/licenses/>.

import sys, os, math
import numpy as np
from Vector import Vector
try :
    import ogr, gdal
except :
    from osgeo import ogr, gdal
from collections import *

[docs]class Segmentation(Vector): """ Vector class inherits the super vector class properties. This class create the final shapefile : Cartography on a input segmentation by decision tree. The output classname are (**out_class_name** variable): - Vegetation non naturelle - Vegetation semi-naturelle - Herbacees - Ligneux - Ligneux mixtes - Ligneux denses - Forte phytomasse - Moyenne phytomasse - Faible phytomasse :param vector_used: Input/Output shapefile to clip (path) :type vector_used: str :param vector_cut: Area shapefile (path) :type vector_cut: str :param output_file: Output shapefile cartography. This path is the same than the segmentation path. :type output_file: str :param out_class_name: List of output class name :type out_class_name: list of str :param out_threshold: List of output threshold :type out_threshold: list of str :param max_...: Biomass and density (IDM, SFS index) maximum :type max_...: float :param class_tab_final: Final decision tree table :type class_tab_final: dict :param stats_rpg_tif: Dictionnary with pxl percent in every polygon of class RPG. Mainly 'Maj_count' (Majority Value) and 'Maj_count_perc' (Majority Percent) :type stats_rpg_tif: dict """ def __init__(self, used, cut): """Create a new 'Clustering' instance """ Vector.__init__(self, used, cut) self.stats_rpg_tif = {} self.output_file = 'Final_classification.shp' self.out_class_name = [] self.out_threshold = [] self.class_tab_final = defaultdict(list) self.max_wood_idm = 0 self.max_wood_sfs = 0 self.max_bio = 0
[docs] def create_cartography(self, out_fieldnames, out_fieldtype): """ Function to create a output shapefile. In this output file, there is the final cartography. With output defined field names and field type in the main process. :param out_fieldnames: List of output field names :type out_fieldnames: list of str :param out_fieldtype: List of outpu field type :type out_fieldtype: list of str """ shp_ogr_ds = self.data_source shp_ogr = self.data_source.GetLayer() # Projection # Import input shapefile projection srsObj = shp_ogr.GetSpatialRef() # Conversion to syntax ESRI srsObj.MorphToESRI() ## Remove the output final shapefile if it exists self.vector_used = self.output_file if os.path.exists(self.vector_used): self.data_source.GetDriver().DeleteDataSource(self.vector_used) out_ds = self.data_source.GetDriver().CreateDataSource(self.vector_used) if out_ds is None: print('Could not create file') sys.exit(1) # Specific output layer out_layer = out_ds.CreateLayer(str(self.vector_used), srsObj, geom_type=ogr.wkbMultiPolygon) # Add new fields # To add RPG_CODE field out_fieldnames.insert(2,'RPG_CODE') out_fieldtype.insert(2,ogr.OFTInteger) # To add the others for i in range(0, len(out_fieldnames)): fieldDefn = ogr.FieldDefn(str(out_fieldnames[i]), out_fieldtype[i]) out_layer.CreateField(fieldDefn) # Add 2 fields to convert class string in code to confusion matrix fieldDefn = ogr.FieldDefn('FBPHY_CODE', ogr.OFTInteger) out_layer.CreateField(fieldDefn) out_fieldnames.append('FBPHY_CODE') fieldDefn = ogr.FieldDefn('FBPHY_SUB', ogr.OFTInteger) out_layer.CreateField(fieldDefn) out_fieldnames.append('FBPHY_SUB') # Feature for the ouput shapefile featureDefn = out_layer.GetLayerDefn() in_feature = shp_ogr.SetNextByIndex(0) # Polygons initialisation in_feature = shp_ogr.GetNextFeature() # Loop on input polygons while in_feature: geom = in_feature.GetGeometryRef() # Extract input geometry # Create a new polygon out_feature = ogr.Feature(featureDefn) # Set the polygon geometry and attribute out_feature.SetGeometry(geom) # Set the existing ID out_feature.SetField(out_fieldnames[0], in_feature.GetField(self.field_names[2])) # Set the area out_feature.SetField(out_fieldnames[1], geom.GetArea()/10000) # 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'] out_feature.SetField('RPG_CODE', int(recouv_crops_RPG)) # Set the others polygons fields with the decision tree dictionnary for i in range(3, len(out_fieldnames)): # If list stopped it on the second level, complete by empty case if len(self.class_tab_final[in_feature.GetFID()]) < len(out_fieldnames)-3 and \ self.class_tab_final[in_feature.GetFID()] != []: self.class_tab_final[in_feature.GetFID()].insert(len(self.class_tab_final[in_feature.GetFID()])-3,'') # To 3rd level self.class_tab_final[in_feature.GetFID()].insert(len(self.class_tab_final[in_feature.GetFID()])-3,0) # To degree try: # To the confusion matrix, replace level ++ by level -- if i == len(out_fieldnames)-1: if self.class_tab_final[in_feature.GetFID()][i-3] == 6: # Crops to artificial vegetation self.class_tab_final[in_feature.GetFID()][i-4] = 0 # if self.class_tab_final[in_feature.GetFID()][i-3] == 2: # # Grassland to natural vegetation # self.class_tab_final[in_feature.GetFID()][i-3] = 1 if self.class_tab_final[in_feature.GetFID()][i-3] > 7: # Phytomass to natural vegetation self.class_tab_final[in_feature.GetFID()][i-4] = 1 out_feature.SetField(str(out_fieldnames[i]), self.class_tab_final[in_feature.GetFID()][i-3]) except: # pass for i in range(3, len(out_fieldnames)-3): out_feature.SetField(str(out_fieldnames[i]), 'Undefined') out_feature.SetField('FBPHY_CODE', 255) out_feature.SetField('FBPHY_SUB', 255) # sys.exit(1) # Append polygon to the output shapefile out_layer.CreateFeature(out_feature) # Destroy polygons out_feature.Destroy() in_feature.Destroy() # Next polygon in_feature = shp_ogr.GetNextFeature() # Close data out_ds.Destroy()
[docs] def decision_tree(self, combin_tree): """ Function to build the decision tree. Taking account output threshold and input class name. :param combin_tree: Decision tree combination :type combin_tree: list of number class name """ # Combination tree on a sentence. Every sentence will be in a table. cond_tab = [] for ct in combin_tree: cond_a = '' # Condition Term c = 0 while c < len(ct): # Loop on tree combinaison if self.out_threshold[ct[c]] =='': # For interval condition cond_a = cond_a + 'self.stats_dict[ind_stats][' + str(ct[c]/2) + ']' +\ self.out_threshold[ct[c]-1].replace('>', '<=') + \ ' and self.stats_dict[ind_stats][' + str(ct[c]/2) + ']' +\ self.out_threshold[ct[c]+1].replace('<', '>=') + ' and ' else: cond_a = cond_a + 'self.stats_dict[ind_stats][' + str(ct[c]/2) + ']' +\ self.out_threshold[ct[c]] + ' and ' c = c + 1 cond_tab.append(cond_a[:-5]) # Remove the last 'and' # Loop on every value for ind_stats in range(len(self.stats_dict)): # Loop on decision tree combination. for cond in cond_tab: # Add class name in the output table try: if eval(cond): self.class_tab_final[ind_stats] = [self.out_class_name[s] \ for s in combin_tree[cond_tab.index(cond)]] + \ [combin_tree[cond_tab.index(cond)][len(combin_tree[cond_tab.index(cond)])-1]] + \ [combin_tree[cond_tab.index(cond)][len(combin_tree[cond_tab.index(cond)])-1]] except NameError: # If there is 'nan' in the table statistics if eval(cond.replace('nan','-10000')):# If there is 'nan' in the table statistics self.class_tab_final[ind_stats] = [self.out_class_name[s] \ for s in combin_tree[cond_tab.index(cond)]] + \ [combin_tree[cond_tab.index(cond)][len(combin_tree[cond_tab.index(cond)])-1]] + \ [combin_tree[cond_tab.index(cond)][len(combin_tree[cond_tab.index(cond)])-1]]
[docs] def compute_biomass_density(self, method='SEATH'): """ Function to compute the biomass and density distribution. It returns threshold of biomass level. :param method: Classification method used. It can set 'SEATH' (by default) or 'RF' :type method: str """ if method == 'SEATH': distri = [v[1:] for k, v in self.stats_dict.items() if eval('v[0]' + self.out_threshold[1])] distri_bio = [] distri_den = [] for b in distri: if eval('b[0]' + self.out_threshold[2]) and b[len(b)-1] != float('inf') and b[len(b)-1] != float('nan') and b[len(b)-1] < 1: distri_bio.append(b) else: distri_den.append(b) elif method == 'RF': distri = [v[1:] for k, v in self.stats_dict.items() if not self.out_threshold[k] in [0,6,7]] distri_bio = [] distri_den = [] for b in distri: if self.out_threshold[distri.index(b)] in [1,2,8,9,10] and b[len(b)-1] != -10000 and b[len(b)-1] < 1: distri_bio.append(b) else: distri_den.append(b) # Set this variable used normally to define threshold of the classification with SEATH method self.out_threshold = [] # Transpose table t_distri_bio = list(map(list, zip(*distri_bio))) t_distri_den = list(map(list, zip(*distri_den))) # Biomass threshold 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)) # 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])
[docs] def append_scale(self, select_class, form): """ Function to complete the 'class_tab_final' list with density and biomass information. This list will be used to build the final shapefile. :param select_class: Class name to add degree :type select_class: str :param form: Formula to add degree :type form: str """ 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: self.class_tab_final[ind_stats].insert(len(self.class_tab_final[ind_stats])-2,eval(form)) # To add phytomasse try : if self.class_tab_final[ind_stats][self.class_tab_final[ind_stats].index(eval(form))-1] == '' and \ self.class_tab_final[ind_stats][self.class_tab_final[ind_stats].index(eval(form))-2] != '': # dict[][A.index(real_pourcent)-1] == '' and dict[][A.index(real_pourcent)-2] != '' # Print phytomasse class in the tab because of self.in_class_name in the Processing class if not eval(form + self.out_threshold[0]) and not eval(form + self.out_threshold[2]): self.class_tab_final[ind_stats][self.class_tab_final[ind_stats].index(eval(form))-1] = self.out_class_name[9] self.class_tab_final[ind_stats][len(self.class_tab_final[ind_stats])-1] = 9 elif eval(form + self.out_threshold[0]): self.class_tab_final[ind_stats][self.class_tab_final[ind_stats].index(eval(form))-1] = self.out_class_name[8] self.class_tab_final[ind_stats][len(self.class_tab_final[ind_stats])-1] = 8 elif eval(form + self.out_threshold[2]): self.class_tab_final[ind_stats][self.class_tab_final[ind_stats].index(eval(form))-1] = self.out_class_name[10] self.class_tab_final[ind_stats][len(self.class_tab_final[ind_stats])-1] = 10 except ValueError: pass except IndexError: pass