Commit a639f083 authored by Commandre Benjamin's avatar Commandre Benjamin
Browse files

Nettoyage du code + Optimisation

1 merge request!1Develop
Showing with 590 additions and 577 deletions
+590 -577
This diff is collapsed.
This diff is collapsed.
...@@ -3,17 +3,17 @@ ...@@ -3,17 +3,17 @@
# #
# This file is part of PHYMOBAT 1.2. # This file is part of PHYMOBAT 1.2.
# Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS) # Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS)
# #
# PHYMOBAT 1.2 is free software: you can redistribute it and/or modify # 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 # it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or # the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version. # (at your option) any later version.
# #
# PHYMOBAT 1.2 is distributed in the hope that it will be useful, # PHYMOBAT 1.2 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of # but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details. # GNU General Public License for more details.
# #
# You should have received a copy of the GNU General Public License # 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/>. # along with PHYMOBAT 1.2. If not, see <http://www.gnu.org/licenses/>.
...@@ -31,9 +31,9 @@ except : ...@@ -31,9 +31,9 @@ except :
class Sample(Vector): class Sample(Vector):
""" """
Vector class inherits the super vector class properties. Vector class inherits the super vector class properties.
This class create training sample. This class create training sample.
@param vector_used: Input/Output shapefile to clip (path) @param vector_used: Input/Output shapefile to clip (path)
@type vector_used: str @type vector_used: str
...@@ -45,32 +45,32 @@ class Sample(Vector): ...@@ -45,32 +45,32 @@ class Sample(Vector):
@param vector_val: Output shapefile to validate the futur classification @param vector_val: Output shapefile to validate the futur classification
@type vector_val: str @type vector_val: str
@opt: Refer to the Vector class @opt: Refer to the Vector class
""" """
def __init__(self, used, cut, nb_sample, **opt): def __init__(self, used, cut, nb_sample, **opt):
""" """
Create a new 'Sample' instance Create a new 'Sample' instance
""" """
self.logger = Outils.Log('Log', 'Sample') self.logger = Outils.Log('Log', 'Sample')
Vector.__init__(self, used, cut, **opt) Vector.__init__(self, used, cut, **opt)
self._nb_sample = nb_sample self._nb_sample = nb_sample
self.vector_val = '' self.vector_val = ''
def create_sample(self, **kwargs): def create_sample(self, **kwargs):
""" """
Function to create a sample shapefile of a specific class Function to create a sample shapefile of a specific class
:kwargs: **fieldname** (list of str) - Fieldname in the input shapefile (if the user want select polygons of the class names specific) :kwargs: **fieldname** (list of str) - Fieldname in the input shapefile (if the user want select polygons of the class names specific)
**class** (list of str) - class names in the input shapefile (with fieldname index). **class** (list of str) - class names in the input shapefile (with fieldname index).
Can use one or several classes like this --> example : [classname1, classname2, ...] Can use one or several classes like this --> example : [classname1, classname2, ...]
""" """
kw_field = kwargs['fieldname'] if kwargs.get('fieldname') else '' kw_field = kwargs['fieldname'] if kwargs.get('fieldname') else ''
kw_classes = kwargs['class'] if kwargs.get('class') else '' kw_classes = kwargs['class'] if kwargs.get('class') else ''
...@@ -78,7 +78,7 @@ class Sample(Vector): ...@@ -78,7 +78,7 @@ class Sample(Vector):
self.logger.debug("kw_classes : {0}".format(kw_classes)) self.logger.debug("kw_classes : {0}".format(kw_classes))
self.output_folder = "{0}/{1}".format(self.output_folder, os.path.basename(os.path.dirname(self.vector_used))) self.output_folder = "{0}/{1}".format(self.output_folder, os.path.basename(os.path.dirname(self.vector_used)))
if not os.path.exists(self.output_folder) : if not os.path.exists(self.output_folder) :
os.mkdir(self.output_folder) os.mkdir(self.output_folder)
...@@ -89,7 +89,7 @@ class Sample(Vector): ...@@ -89,7 +89,7 @@ class Sample(Vector):
else: else:
# The random sample without class name selected # The random sample without class name selected
random_sample = np.array(random.sample(range(self.data_source.GetLayer().GetFeatureCount()), self._nb_sample)) random_sample = np.array(random.sample(range(self.data_source.GetLayer().GetFeatureCount()), self._nb_sample))
self.logger.debug("random_sample : {0}".format(random_sample)) self.logger.debug("random_sample : {0}".format(random_sample))
# Output shapefile of the sample's polygons (path) # Output shapefile of the sample's polygons (path)
...@@ -102,61 +102,51 @@ class Sample(Vector): ...@@ -102,61 +102,51 @@ class Sample(Vector):
# Fill and create the sample shapefile # Fill and create the sample shapefile
self.fill_sample(self.vector_used, random_sample[:round(len(random_sample)/2)]) self.fill_sample(self.vector_used, random_sample[:round(len(random_sample)/2)])
# Output shapefile of the validate polygon (path) # Output shapefile of the validate polygon (path)
self.vector_val = "{0}val.shp".format(self.vector_used[:-6]) self.vector_val = "{0}val.shp".format(self.vector_used[:-6])
# Fill and create the validate polygons shapefile # Fill and create the validate polygons shapefile
self.fill_sample(self.vector_val, random_sample[round(len(random_sample)/2):]) self.fill_sample(self.vector_val, random_sample[round(len(random_sample)/2):])
def select_random_sample(self, kw_field, kw_classes): def select_random_sample(self, kw_field, kw_classes):
""" """
Function to select id with class name specific only. Function to select id with class name specific only.
This function is used in :func:`create_sample` This function is used in :func:`create_sample`
@param kw_field: Field name in the input shapefile @param kw_field: Field name in the input shapefile
@type kw_field: str @type kw_field: str
@param kw_classes: Class names in the input shapefile like this --> 'classname1, classname2' @param kw_classes: Class names in the input shapefile like this --> 'classname1, classname2'
@type kw_classes: str @type kw_classes: str
@returns: list -- variable **select_id**, List of id with a class name specific. @returns: list -- variable **select_id**, List of id with a class name specific.
""" """
# Convert string in a list. For that, it remove # Convert string in a list. For that, it remove
# space and clip this string with comma (Add everywhere if the script modified # space and clip this string with comma (Add everywhere if the script modified
# because the process work with a input string chain) # because the process work with a input string chain)
kw_classes = kw_classes.replace(' ','').split(',') kw_classes = kw_classes.replace(' ','').split(',')
self.logger.debug('kw_classes : {0}'.format(kw_classes))
# List of class name id # List of class name id
select_id = [] select_id = []
shp_ogr = self.data_source.GetLayer() layer = self.data_source.GetLayer()
# Loop on input polygons assert layer.GetFeatureCount() != 0, "Le shapefile : {0} est vide.".format(self.data_source.GetLayer().GetName())
in_feature = shp_ogr.SetNextByIndex(0) # Initialization
in_feature = shp_ogr.GetNextFeature() for i in range(layer.GetFeatureCount()) :
feature = layer.GetFeature(i)
while in_feature:
# if polygon is a defined class name if feature.GetFieldIndex(self.field_names[self.field_names.index(kw_field)]) != -1 \
## .replace('0','') to remove '0' in front of for example '1' (RPG -> '01') and not feature.GetField(self.field_names[self.field_names.index(kw_field)]) is None :
table_name_class = in_feature.GetField(self.field_names[self.field_names.index(kw_field)]) if feature.GetField(self.field_names[self.field_names.index(kw_field)]).replace('0','') in kw_classes :
# To avoid that the process crashed this part of the algorithm will be launch select_id.append(i)
# if the field is contains characters
if table_name_class != None : assert len(select_id) > 0, "Aucun polygone ne correspond."
if in_feature.GetField(self.field_names[self.field_names.index(kw_field)]).replace('0','') in kw_classes:
# Add id in the extract list
select_id.append(in_feature.GetFID())
in_feature.Destroy()
in_feature = shp_ogr.GetNextFeature()
return select_id return select_id
def fill_sample(self, output_sample, polygon, **opt): def fill_sample(self, output_sample, polygon, **opt):
""" """
Function to fill and create the output sample shapefile. Function to fill and create the output sample shapefile.
...@@ -165,15 +155,15 @@ class Sample(Vector): ...@@ -165,15 +155,15 @@ class Sample(Vector):
@param output_sample: Path of the output shapefile @param output_sample: Path of the output shapefile
@type output_sample: str @type output_sample: str
@param polygon: Identity of the selected random polygons. @param polygon: Identity of the selected random polygons.
If this variable = 0, the processing will take all polygons If this variable = 0, the processing will take all polygons
@type polygon: list or int @type polygon: list or int
@opt: **add_fieldname** (int) - Variable to kown if add a field. By default non (0), if it have to add (1) @opt: **add_fieldname** (int) - Variable to kown if add a field. By default non (0), if it have to add (1)
**fieldname** (str) - Fieldname to add in the input shapefile **fieldname** (str) - Fieldname to add in the input shapefile
**class** (int) - class names in integer to add in the input shapefile **class** (int) - class names in integer to add in the input shapefile
""" """
...@@ -181,56 +171,56 @@ class Sample(Vector): ...@@ -181,56 +171,56 @@ class Sample(Vector):
add_field = opt['add_fieldname'] if opt.get('add_fieldname') else 0 add_field = opt['add_fieldname'] if opt.get('add_fieldname') else 0
opt_field = opt['fieldname'] if opt.get('fieldname') else '' opt_field = opt['fieldname'] if opt.get('fieldname') else ''
opt_class = opt['class'] if opt.get('class') else 0 opt_class = opt['class'] if opt.get('class') else 0
shp_ogr = self.data_source.GetLayer() shp_ogr = self.data_source.GetLayer()
# To take all polygon # To take all polygon
if type(polygon) == int: if type(polygon) == int:
polygon = range(shp_ogr.GetFeatureCount()) polygon = range(shp_ogr.GetFeatureCount())
# Projection # Projection
# Import input shapefile projection # Import input shapefile projection
srsObj = shp_ogr.GetSpatialRef() srsObj = shp_ogr.GetSpatialRef()
# Conversion to syntax ESRI # Conversion to syntax ESRI
srsObj.MorphToESRI() srsObj.MorphToESRI()
## Remove the output shapefile if it exists ## Remove the output shapefile if it exists
if os.path.exists(output_sample): if os.path.exists(output_sample):
self.data_source.GetDriver().DeleteDataSource(output_sample) self.data_source.GetDriver().DeleteDataSource(output_sample)
out_ds = self.data_source.GetDriver().CreateDataSource(output_sample) out_ds = self.data_source.GetDriver().CreateDataSource(output_sample)
if out_ds is None: if out_ds is None:
self.logger.error('Could not create file') self.logger.error('Could not create file')
sys.exit(1) sys.exit(1)
# Specific output layer # Specific output layer
out_layer = out_ds.CreateLayer(output_sample, srsObj, geom_type=ogr.wkbMultiPolygon) out_layer = out_ds.CreateLayer(output_sample, srsObj, geom_type=ogr.wkbMultiPolygon)
# Add existing fields # Add existing fields
for i in range(0, len(self.field_names)): for i in range(0, len(self.field_names)):
# use the input FieldDefn to add a field to the output # use the input FieldDefn to add a field to the output
fieldDefn = shp_ogr.GetFeature(0).GetFieldDefnRef(self.field_names[i]) fieldDefn = shp_ogr.GetFeature(0).GetFieldDefnRef(self.field_names[i])
out_layer.CreateField(fieldDefn) out_layer.CreateField(fieldDefn)
# In Option : Add a integer field # In Option : Add a integer field
if add_field == 1: if add_field == 1:
new_field = ogr.FieldDefn(opt_field, 0) new_field = ogr.FieldDefn(opt_field, 0)
out_layer.CreateField(new_field) out_layer.CreateField(new_field)
# Feature for the ouput shapefile # Feature for the ouput shapefile
featureDefn = out_layer.GetLayerDefn() featureDefn = out_layer.GetLayerDefn()
# Loop on the input elements # Loop on the input elements
# Create a existing polygons in random list # Create a existing polygons in random list
for cnt in polygon: for cnt in polygon:
# Select input polygon by id # Select input polygon by id
in_feature = shp_ogr.SetNextByIndex(cnt) in_feature = shp_ogr.SetNextByIndex(cnt)
in_feature = shp_ogr.GetNextFeature() in_feature = shp_ogr.GetNextFeature()
# Extract input geometry # Extract input geometry
geom = in_feature.GetGeometryRef() geom = in_feature.GetGeometryRef()
# Create a new polygon # Create a new polygon
out_feature = ogr.Feature(featureDefn) out_feature = ogr.Feature(featureDefn)
...@@ -249,10 +239,10 @@ class Sample(Vector): ...@@ -249,10 +239,10 @@ class Sample(Vector):
# Append polygon to the output shapefile # Append polygon to the output shapefile
out_layer.CreateFeature(out_feature) out_layer.CreateFeature(out_feature)
# Destroy polygons # Destroy polygons
out_feature.Destroy() out_feature.Destroy()
in_feature.Destroy() in_feature.Destroy()
# Close data # Close data
out_ds.Destroy() out_ds.Destroy()
\ No newline at end of file
...@@ -3,39 +3,34 @@ ...@@ -3,39 +3,34 @@
# #
# This file is part of PHYMOBAT 1.2. # This file is part of PHYMOBAT 1.2.
# Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS) # Copyright 2016 Sylvio Laventure (IRSTEA - UMR TETIS)
# #
# PHYMOBAT 1.2 is free software: you can redistribute it and/or modify # 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 # it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or # the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version. # (at your option) any later version.
# #
# PHYMOBAT 1.2 is distributed in the hope that it will be useful, # PHYMOBAT 1.2 is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of # but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details. # GNU General Public License for more details.
# #
# You should have received a copy of the GNU General Public License # 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/>. # along with PHYMOBAT 1.2. If not, see <http://www.gnu.org/licenses/>.
import os, sys import os, sys
import subprocess import subprocess
import numpy as np import numpy as np
import osgeo from osgeo import ogr
# try :
# import ogr, gdal
# except :
# from osgeo import ogr, gdal
import otbApplication as otb import otbApplication as otb
from collections import * from collections import *
import Outils import Outils
from RasterSat_by_date import RasterSat_by_date from RasterSat_by_date import RasterSat_by_date
class Vector(): class Vector():
""" """
Vector class to extract a area, vector data and zonal statistic Vector class to extract a area, vector data and zonal statistic
@param vector_used: Input/Output shapefile to clip (path) @param vector_used: Input/Output shapefile to clip (path)
@type vector_used: str @type vector_used: str
...@@ -59,10 +54,10 @@ class Vector(): ...@@ -59,10 +54,10 @@ class Vector():
@opt: **Remove** (int) - For the remove_shp variable @opt: **Remove** (int) - For the remove_shp variable
""" """
def __init__(self, used, cut, **opt): def __init__(self, used, cut, **opt):
""" """
Create a new 'Vector' instance Create a new 'Vector' instance
""" """
if not hasattr(self, "logger") : if not hasattr(self, "logger") :
...@@ -74,19 +69,19 @@ class Vector(): ...@@ -74,19 +69,19 @@ class Vector():
self.vector_cut = cut self.vector_cut = cut
self.vector_used = used self.vector_used = used
self.remove_shp = opt['Remove'] if opt.get('Remove') else 0 self.remove_shp = opt['Remove'] if opt.get('Remove') else 0
self.stats_dict = defaultdict(list) self.stats_dict = defaultdict(list)
def clip_vector(self, output_folder): def clip_vector(self, output_folder):
""" """
Function to clip a vector with a vector Function to clip a vector with a vector
""" """
if not os.path.exists(output_folder) : if not os.path.exists(output_folder) :
os.makedirs(output_folder) os.makedirs(output_folder)
outclip = "{0}/Clip_{1}".format(output_folder, self.vector_name) outclip = "{0}/Clip_{1}".format(output_folder, self.vector_name)
if not os.path.exists(outclip) or self.remove_shp == 1: if not os.path.exists(outclip) or self.remove_shp == 1:
self.logger.info('Clip of {0}'.format(self.vector_name)) self.logger.info('Clip of {0}'.format(self.vector_name))
...@@ -100,43 +95,43 @@ class Vector(): ...@@ -100,43 +95,43 @@ class Vector():
# Replace input filename by output filename # Replace input filename by output filename
self.vector_used = outclip self.vector_used = outclip
def vector_data(self): def vector_data(self):
""" """
Function to extract vector layer information Function to extract vector layer information
""" """
try: try:
self.data_source = osgeo.ogr.GetDriverByName('ESRI Shapefile').Open(self.vector_used, 0) self.data_source = ogr.Open(self.vector_used)
self.logger.info('Shapefile opening : {0}'.format(self.data_source.GetLayer().GetName())) self.logger.info('Shapefile opening : {0}'.format(self.data_source.GetLayer().GetName()))
except : except :
self.logger.error('Could not open file') self.logger.error('Could not open file')
sys.exit(1) sys.exit(1)
# List of field name # List of field name
self.field_names = [self.data_source.GetLayer().GetLayerDefn().GetFieldDefn(l).GetName() \ self.field_names = [self.data_source.GetLayer().GetLayerDefn().GetFieldDefn(l).GetName() \
for l in range(self.data_source.GetLayer().GetLayerDefn().GetFieldCount())] for l in range(self.data_source.GetLayer().GetLayerDefn().GetFieldCount())]
def close_data(self): def close_data(self):
""" """
Function to remove allocate memory Function to remove allocate memory
""" """
self.logger.info('Shapefile closing : {0}'.format(self.data_source.GetLayer().GetName()))
# Close data sources # Close data sources
self.data_source.Destroy() self.data_source.Destroy()
self.logger.info('Shapefile closing : {0}'.format(self.data_source.GetLayer().GetName()))
def zonal_stats(self, liste_chemin): def zonal_stats(self, liste_chemin):
""" """
Function to compute the mean in every polygons on a list images with otb Function to compute the mean in every polygons on a list images with otb
:param liste_chemin : List input image path @param liste_chemin : List input image path
:type liste_chemin: list(string) @type liste_chemin : list(string)
""" """
self.logger.info("Compute stats 'mean' on {0}".format(os.path.split(self.vector_used)[1])) self.logger.info("Compute stats 'mean' on {0}".format(os.path.split(self.vector_used)[1]))
zonal_stats = otb.Registry.CreateApplication("ZonalStatistics") zonal_stats = otb.Registry.CreateApplication("ZonalStatistics")
zonal_stats.SetParameterStringList("il", liste_chemin) zonal_stats.SetParameterStringList("il", liste_chemin)
zonal_stats.SetParameterStringList("vec", [self.vector_used]) zonal_stats.SetParameterStringList("vec", [self.vector_used])
...@@ -159,15 +154,15 @@ class Vector(): ...@@ -159,15 +154,15 @@ class Vector():
for idx, r in enumerate(resultats) : for idx, r in enumerate(resultats) :
self.stats_dict[idx] = r.tolist() self.stats_dict[idx] = r.tolist()
self.logger.info("End of stats 'mean' on {0}".format(os.path.split(self.vector_used)[1])) self.logger.info("End of stats 'mean' on {0}".format(os.path.split(self.vector_used)[1]))
def zonal_stats_pp(self, inraster): def zonal_stats_pp(self, inraster):
""" """
A zonal statistics ++ to dertermine pxl percent in every polygon A zonal statistics ++ to dertermine pxl percent in every polygon
:param inraster: Input image path :param inraster: Input image path
:type inraster: str :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) :returns: dict -- **p_stats** : dictionnary with pxl percent in every polygon. Mainly 'Maj_count' (Majority Value) and 'Maj_count_perc' (Majority Percent)
""" """
...@@ -191,30 +186,30 @@ class Vector(): ...@@ -191,30 +186,30 @@ class Vector():
dico['Maj_count'] = float(liste_item[1]) dico['Maj_count'] = float(liste_item[1])
dico['Maj_count_perc'] = float(liste_item[2]) dico['Maj_count_perc'] = float(liste_item[2])
p_stats.append(dico) p_stats.append(dico)
return p_stats return p_stats
def layer_rasterization(self, raster_head, attribute_r, **kwargs): def layer_rasterization(self, raster_head, attribute_r, **kwargs):
""" """
Function to rasterize a vector using OTB. Function to rasterize a vector using OTB.
@param raster_head: Raster path that will look like the final raster of the rasterization @param raster_head: Raster path that will look like the final raster of the rasterization
@type raster_head: str @type raster_head: str
@param attribute_r: Field name of the shapefile that contains class names @param attribute_r: Field name of the shapefile that contains class names
@type attribute_r: str @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... @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 @returns: str -- **valid_raster** : output raster path from the rasterization
""" """
valid_raster = "{0}.TIF".format(self.vector_used[:-4])# Name of the output raster valid_raster = "{0}.TIF".format(self.vector_used[:-4])# Name of the output raster
self.logger.debug("valid_raster : {0}".format(valid_raster)) self.logger.debug("valid_raster : {0}".format(valid_raster))
if os.path.exists(valid_raster): if os.path.exists(valid_raster):
os.remove(valid_raster) os.remove(valid_raster)
layerRasterization = otb.Registry.CreateApplication("Rasterization") layerRasterization = otb.Registry.CreateApplication("Rasterization")
layerRasterization.SetParameterString("in", self.vector_used) layerRasterization.SetParameterString("in", self.vector_used)
layerRasterization.SetParameterString("out", valid_raster) layerRasterization.SetParameterString("out", valid_raster)
...@@ -223,5 +218,5 @@ class Vector(): ...@@ -223,5 +218,5 @@ class Vector():
layerRasterization.SetParameterString("mode.attribute.field", attribute_r) layerRasterization.SetParameterString("mode.attribute.field", attribute_r)
layerRasterization.ExecuteAndWriteOutput() layerRasterization.ExecuteAndWriteOutput()
return valid_raster return valid_raster
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment