Commit 5c9bdcd67ebcb62e5ccdbdbe26eb4f5601265870

Authored by Louis Baetens
1 parent b1158dcb
Exists in master

RMV unused directories

ALCD/parameters_files/global_parameters.json
... ... @@ -103,11 +103,11 @@
103 103 "training_proportion": "0.9"
104 104 },
105 105 "user_choices": {
106   - "clear_date": "20180823",
107   - "current_date": "20180813",
108   - "location": "Alta_Floresta_Brazil",
109   - "main_dir": "/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/Alta_Floresta_Brazil_21LWK_20180813",
110   - "raw_img": "Alta_Floresta_Brazil_bands.tif",
111   - "tile": "21LWK"
  106 + "clear_date": "20160427",
  107 + "current_date": "20160417",
  108 + "location": "Marrakech",
  109 + "main_dir": "/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/Marrakech_29RPQ_20160417",
  110 + "raw_img": "Marrakech_bands.tif",
  111 + "tile": "29RPQ"
112 112 }
113 113 }
114 114 \ No newline at end of file
... ...
ALCD/unused/GDAL_extract_part.py
... ... @@ -1,41 +0,0 @@
1   -import os
2   -import gdal
3   -
4   -
5   -def extract_part(in_file, out_file, left_side = 0.3,right_side = 0.4,top_side = 0.7,bottom_side = 0.9):
6   - ''' Extract a part of a tif, given the positions of the rectangle sides
7   - given in percentage of total image
8   - '''
9   - #~ input_file = '/mnt/data/home/baetensl/OTB_codes/OTB_in/B04_B03_B02_B10.tif'
10   - #~ output_file = '/mnt/data/home/baetensl/OTB_codes/OTB_in/cutted2.tif'
11   -
12   - # get the corners of the original image
13   - src = gdal.Open(input_file)
14   - ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform() #upper left corner x and y
15   - lrx = ulx + (src.RasterXSize * xres)
16   - lry = uly + (src.RasterYSize * yres) #low right corner x and y
17   -
18   - x_width = lrx-ulx
19   - y_height = lry-uly
20   -
21   - #~ left_side = 0.3
22   - #~ right_side = 0.4
23   - #~ top_side = 0.7
24   - #~ bottom_side = 0.9
25   -
26   - cut_ulx = ulx+left_side*x_width
27   - cut_lrx = ulx+right_side*x_width
28   -
29   - cut_uly = uly+top_side*y_height
30   - cut_lry = uly+bottom_side*y_height
31   -
32   - print([cut_ulx, cut_lrx, cut_uly, cut_lry])
33   -
34   - command = 'gdal_translate -projwin {} {} {} {} {} {} -of GTiff'.format(cut_ulx, cut_uly, cut_lrx, cut_lry, input_file, output_file)
35   - os.system(command)
36   -
37   -if __name__ == '__main__':
38   - input_file = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/In_data/Image/all_bands.tif'
39   - output_file = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/In_data/Image/cutted.tif'
40   - extract_part(in_file = input_file, out_file = output_file, left_side = 0.3,right_side = 0.4,top_side = 0.7,bottom_side = 0.9)
41   -
ALCD/unused/OTB_dimension_reduction.py
... ... @@ -1,48 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -import os
4   -import os.path as op
5   -import json
6   -import otbApplication
7   -import csv
8   -import subprocess
9   -import xml.etree.ElementTree as ET
10   -
11   -
12   -def reduction(in_tif, out_tif):
13   - # The following line creates an instance of the DimensionalityReduction application
14   - DimensionalityReduction = otbApplication.Registry.CreateApplication("DimensionalityReduction")
15   -
16   - # The following lines set all the application parameters:
17   - DimensionalityReduction.SetParameterString("in", str(in_tif))
18   -
19   - DimensionalityReduction.SetParameterString("out", str(out_tif))
20   -
21   - DimensionalityReduction.SetParameterString("method","pca")
22   - DimensionalityReduction.SetParameterString("normalize","True")
23   - DimensionalityReduction.SetParameterInt("nbcomp",3)
24   - DimensionalityReduction.SetParameterString("rescale.outmin",'-1.0')
25   - DimensionalityReduction.SetParameterString("rescale.outmax",'1.0')
26   -
27   - DimensionalityReduction.UpdateParameters()
28   - # The following line execute the application
29   - DimensionalityReduction.ExecuteAndWriteOutput()
30   -
31   -
32   -
33   -def main():
34   - main_dir = '/mnt/data/home/baetensl/classification_clouds/Data/Arles_20171002/In_data/Image'
35   -
36   - in_tif = op.join(main_dir, 'Arles_bands.tif')
37   - out_tif = op.join(main_dir, 'dimension_reduced.tif')
38   - reduction(in_tif, out_tif)
39   -
40   -
41   -
42   -if __name__ == '__main__':
43   - main()
44   -
45   -
46   -
47   -
48   -
ALCD/unused/OTB_resample.py
... ... @@ -1,128 +0,0 @@
1   -#!/usr/bin/python
2   -# -*- coding: utf-8 -*-
3   -
4   -import sys
5   -import os
6   -import os.path as op
7   -import otbApplication
8   -import subprocess
9   -
10   -
11   -def morpho_operation(in_tif, out_tif, morph_type, radius, foreground_value = 1.0, background_value = 2.0):
12   - ''' Performs an erosion or a dilatation on the in_tif
13   - Structural element size is radius
14   - '''
15   - print(" Morphological operation")
16   - BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
17   -
18   - BinaryMorphologicalOperation.SetParameterString("in", in_tif)
19   - BinaryMorphologicalOperation.SetParameterString("out", out_tif)
20   - BinaryMorphologicalOperation.SetParameterInt("channel", 1)
21   -
22   - BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", radius)
23   - BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", radius)
24   -
25   - if morph_type == 'erosion':
26   - BinaryMorphologicalOperation.SetParameterString("filter","erode")
27   - BinaryMorphologicalOperation.SetParameterFloat("filter.erode.foreval", foreground_value)
28   - BinaryMorphologicalOperation.SetParameterFloat("filter.erode.backval", background_value)
29   -
30   - elif morph_type == 'dilatation':
31   - BinaryMorphologicalOperation.SetParameterString("filter","dilate")
32   - BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", foreground_value)
33   - BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.backval", background_value)
34   -
35   - BinaryMorphologicalOperation.UpdateParameters()
36   - BinaryMorphologicalOperation.ExecuteAndWriteOutput()
37   - print("Done")
38   -
39   -
40   -def images_intersection(img_tif1, img_tif2, out_img):
41   - ''' Intersection of the two img, output will have the same
42   - value if the original intersect, and 0 elsewhere
43   - '''
44   - print(" Mathematical intersection")
45   - BandMath = otbApplication.Registry.CreateApplication("BandMath")
46   - BandMath.SetParameterStringList("il", [img_tif1, img_tif2])
47   - BandMath.SetParameterString("out", out_img)
48   -
49   - BandMath.SetParameterString("exp", "(im1b1 == im2b1 ? im1b1:0)")
50   -
51   - BandMath.UpdateParameters()
52   - BandMath.ExecuteAndWriteOutput()
53   - print("Done")
54   -
55   -
56   -
57   -def polygonize_tiff(in_tif, mask_tif, out_shp):
58   - ''' Polygonize a tiff to a shapefile
59   - '''
60   - print(" Polygonization")
61   -
62   - out_layer = out_shp.split(op.sep)[-1].replace('.shp', '')
63   -
64   - command = 'gdal_polygonize.py {} -mask {} -f "ESRI Shapefile" {} {} class'.format(in_tif, mask_tif, out_shp, out_layer)
65   - subprocess.call(command, shell=True)
66   -
67   - print("Done")
68   -
69   -def rasterize_shapefile(in_shp, out_tif):
70   - print(" Rasterization")
71   - layer_name = 'layerTest'
72   -
73   - command = 'gdal_rasterize -a class -l points_clouds {in_shp} {out_tif}'.format(layer_name=layer_name, in_shp=in_shp, out_tif = out_tif)
74   - print(command)
75   - subprocess.call(command, shell=True)
76   -
77   - print("Done")
78   -
79   -
80   -
81   -def simplify_geometry(in_shp, out_shp, tolerance = 100):
82   - ''' Simplification of a shapefile
83   - This allows to have lighter polygons, enhancing the rapidty
84   - '''
85   - print(" Geometry simplification")
86   - command = "ogr2ogr {} {} -simplify {}".format(out_shp, in_shp, str(tolerance))
87   - subprocess.call(command, shell=True)
88   - print("Done")
89   -
90   -
91   -def resample(labeled_img, out_shp):
92   - ''' Main function of this file
93   - '''
94   - intermediate_dir = op.join(op.dirname(labeled_img), '..', 'Intermediate')
95   -
96   - dilated_img = op.join(intermediate_dir, 'dilated.tif')
97   - eroded_img = op.join(intermediate_dir, 'eroded.tif')
98   - combined_img = op.join(intermediate_dir, 'combined.tif')
99   - complex_shp = op.join(intermediate_dir, 'complex_shp.shp')
100   -
101   - for filen in [dilated_img, eroded_img, combined_img, complex_shp]:
102   - if op.exists(filen):
103   - os.remove(filen)
104   -
105   - # dilatation direction : to the inside of the clouds, erosion direction : to the inside of background
106   - morpho_operation(labeled_img, dilated_img, morph_type = "dilatation", radius = 4)
107   - morpho_operation(labeled_img, eroded_img, morph_type = "erosion", radius = 10)
108   -
109   - images_intersection(eroded_img, dilated_img, combined_img)
110   -
111   - polygonize_tiff(eroded_img, combined_img, complex_shp)
112   -
113   - simplify_geometry(complex_shp, out_shp, tolerance = 100)
114   -
115   -
116   -def main():
117   - labeled_img = 'labeled_image.tif'
118   - out_shp = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Python_dir/In_data/Masks/output_shp.shp'
119   -
120   - #~ resample(labeled_img, out_shp)
121   -
122   - in_shp = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/In_data/Masks/points_clouds.shp'
123   - out_tif = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/In_data/Masks/points_clouds_raster.tif'
124   -
125   - rasterize_shapefile(in_shp, out_tif)
126   -
127   -if __name__=='__main__':
128   - main()
ALCD/unused/OTB_workflow_save.py
... ... @@ -1,421 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -import os, optparse
4   -import sys
5   -import os.path as op
6   -import json
7   -import glob
8   -import otbApplication
9   -import csv
10   -import pandas as pd
11   -import sqlite3
12   -import subprocess
13   -import xml.etree.ElementTree as ET
14   -from OTB_resample import resample
15   -
16   -
17   -#-------- 0. DIRECTORIES CREATION---------------------
18   -def create_directories(main_dir):
19   - '''
20   - 0. Create the directories for the code to work with
21   - '''
22   - directories = ['', 'In_data', op.join('In_data', 'Masks'), op.join('In_data', 'Image'),
23   - 'Statistics', 'Samples', 'Models', 'Out', 'Other', 'Intermediate']
24   - for sub_dir in directories:
25   - current_dir = op.join(main_dir, sub_dir)
26   - if not os.path.exists(current_dir):
27   - os.makedirs(current_dir)
28   - print(current_dir + ' created')
29   -
30   -
31   -
32   -#-------- 1. PREPROCESSING---------------------
33   -#------------ image statistics
34   -def compute_image_stats(main_dir, raw_img, proceed=True):
35   - '''
36   - 4. Compute the images stats
37   - '''
38   - #otbcli_ComputeImagesStatistics -il in/cutted.tif -out out/images_statistics.xml
39   -
40   - img_stats = op.join(main_dir, 'Statistics', 'img_stats.xml')
41   -
42   -
43   - print " Compute Images Statistics"
44   - ComputeImagesStatistics = otbApplication.Registry.CreateApplication("ComputeImagesStatistics")
45   - ComputeImagesStatistics.SetParameterStringList("il", [ raw_img ])
46   - ComputeImagesStatistics.SetParameterString("out", img_stats)
47   - ComputeImagesStatistics.ExecuteAndWriteOutput()
48   -
49   - print('Done')
50   - return img_stats
51   -
52   -
53   -#----------- samples preprocessing
54   -def compute_samples_stats(main_dir, raw_img, training_shp, proceed=True):
55   - '''
56   - 1. Compute the samples stats
57   - '''
58   - # otbcli_PolygonClassStatistics -in in/cutted.tif -vec in/training_clouds.shp -field class -out out/classes.xml
59   - print " Polygon Classes Statistics"
60   -
61   - class_stats = op.join(main_dir, 'Statistics', 'classes.xml')
62   -
63   - PolygonClassStatistics = otbApplication.Registry_CreateApplication("PolygonClassStatistics")
64   - PolygonClassStatistics.SetParameterString("in", raw_img)
65   - PolygonClassStatistics.SetParameterString("vec", training_shp)
66   - PolygonClassStatistics.SetParameterString("out", class_stats)
67   - PolygonClassStatistics.UpdateParameters()
68   - PolygonClassStatistics.SetParameterStringList("field", ["class"])
69   - PolygonClassStatistics.ExecuteAndWriteOutput()
70   - print('Done')
71   -
72   - return class_stats
73   -
74   -def get_bands_qty(img_stats):
75   - '''
76   - Parse the XML file returned by compute_image_stats to get the
77   - correct number of bands used
78   - This is used by function extract_samples()
79   - '''
80   - tree = ET.parse(img_stats)
81   - root = tree.getroot()
82   -
83   - bands_qty = 0
84   - for k in range(0,len(root)):
85   - if root[k].attrib["name"] == "mean":
86   - bands_qty = len(root[k])
87   - print('There are {} bands used'.format(bands_qty))
88   -
89   - return bands_qty
90   -
91   -
92   -def get_samples_nb(class_stats):
93   - '''
94   - Parse the XML file returned by compute_samples_stats to get the
95   - samples number per class
96   - This is used by function select_samples()
97   - '''
98   - tree = ET.parse(class_stats)
99   - root = tree.getroot()
100   -
101   - nbSamples = []
102   - for k in range(0,len(root)):
103   - if root[k].attrib["name"] == "samplesPerClass":
104   - for c in range(0, len(root[k])):
105   - nbSamples.append(int(root[k][c].attrib["value"]))
106   -
107   - return nbSamples
108   -
109   -def select_samples(main_dir, raw_img, training_shp, class_stats, strategy = "smallest", proceed=True):
110   - '''
111   - 2. Select the samples
112   - '''
113   - #otbcli_SampleSelection -in in/cutted.tif -vec in/training_clouds.shp -instats out/classes.xml -field class -strategy smallest -outrates out/rates.csv -out out/samples.sqlite
114   -
115   - training_samples_location = op.join(main_dir, 'Samples', 'training_samples_location.sqlite')
116   - rates = op.join(main_dir, 'Statistics', 'rates.csv')
117   -
118   - print " Training Samples Selection"
119   - SampleSelection = otbApplication.Registry_CreateApplication("SampleSelection")
120   - SampleSelection.SetParameterString("in", raw_img)
121   - SampleSelection.SetParameterString("vec", training_shp)
122   - SampleSelection.SetParameterString("instats", class_stats)
123   - SampleSelection.SetParameterString("out", training_samples_location)
124   - SampleSelection.SetParameterString("outrates", rates)
125   - SampleSelection.SetParameterString("sampler", "random") #default is periodic
126   -
127   - if strategy == "smallest":
128   - SampleSelection.SetParameterString("strategy", "smallest")
129   - elif strategy == "constant":
130   - SampleSelection.SetParameterString("strategy", "constant")
131   - minSamplesNb = min(get_samples_nb(class_stats))
132   - SampleSelection.SetParameterString("strategy.constant.nb", str(minSamplesNb)) #default is periodic
133   - elif strategy.split('_')[0] == "constant":
134   - SampleSelection.SetParameterString("strategy", "constant")
135   - SampleSelection.SetParameterString("strategy.constant.nb", str(strategy.split('_')[1])) #default is periodic
136   -
137   - SampleSelection.UpdateParameters()
138   - SampleSelection.SetParameterStringList("field", ["class"])
139   - SampleSelection.Execute()
140   - print('Done')
141   -
142   - return training_samples_location
143   -
144   -
145   -
146   -def extract_samples(main_dir, raw_img, training_samples_location, proceed=True):
147   - '''
148   - 3. Extract the samples
149   - '''
150   - #otbcli_SampleExtraction -in in/cutted.tif -vec out/samples.sqlite -outfield prefix -outfield.prefix.name band_ -field class -out out/samples_updates.sqlite
151   -
152   - training_samples_extracted = op.join(main_dir, 'Samples', 'training_samples_extracted.sqlite')
153   -
154   - print " Training Samples Extraction"
155   - SampleExtraction = otbApplication.Registry_CreateApplication("SampleExtraction")
156   - SampleExtraction.SetParameterString("in", raw_img)
157   - SampleExtraction.SetParameterString("vec", training_samples_location)
158   - SampleExtraction.SetParameterString("outfield", "prefix")
159   - SampleExtraction.SetParameterString("outfield.prefix.name", "band_")
160   - SampleExtraction.SetParameterString("out", training_samples_extracted)
161   - SampleExtraction.UpdateParameters()
162   - SampleExtraction.SetParameterStringList("field", [ "class" ])
163   - SampleExtraction.ExecuteAndWriteOutput()
164   - print('Done')
165   -
166   - return training_samples_extracted
167   -
168   -
169   -
170   -def samples_preprocessing(main_dir, raw_img, training_shp, proceed=True):
171   - class_stats = compute_samples_stats(main_dir, raw_img, training_shp, proceed=proceed)
172   - training_samples_location = select_samples(main_dir, raw_img, training_shp, class_stats, strategy="constant_8000", proceed=proceed)
173   - training_samples_extracted = extract_samples(main_dir, raw_img, training_samples_location, proceed=proceed)
174   - #~ split_samples_set()
175   -
176   - return training_samples_extracted
177   -
178   -
179   -#------------ total preprocessing
180   -def total_preprocessing(main_dir, raw_img, training_shp, proceed=True):
181   - img_stats = compute_image_stats(main_dir, raw_img)
182   - training_samples_extracted = samples_preprocessing(main_dir, raw_img, training_shp)
183   - return training_samples_extracted, img_stats
184   -
185   -#-------- 2. MODEL TRAINING ---------------------
186   -
187   -def train_model(main_dir, method, training_samples_extracted, img_stats, shell=True, proceed=True):
188   - '''
189   - 5. Train the model
190   - '''
191   - #otbcli_TrainVectorClassifier -io.vd out/samples_updates.sqlite -cfield class -io.out out/model.rf -classifier rf -feat band_0 band_1 band_2 band_3
192   - #mettre l'option en ligne de command ou en python pour voir le progress
193   -
194   - model_out = op.join(main_dir, 'Models/model.' + method)
195   -
196   - features = ''
197   -
198   - nb_feat = get_bands_qty(img_stats)
199   -
200   - for b in range(nb_feat) :
201   - features += (" band_" + str(b))
202   - features = features[1:]
203   - print(str(features))
204   -
205   -
206   - if proceed == True:
207   - print(" Train Vector Classifier")
208   -
209   - if shell == True:
210   - command = 'otbcli_TrainVectorClassifier -io.vd {} -cfield {} -io.out {} -classifier {} -feat {}'.format(training_samples_extracted, "class", model_out, method, str(features))
211   - print(command)
212   - with open(op.join(main_dir, 'Models', 'model_parameters.json'), 'r') as file:
213   - dico = json.load(file)
214   -
215   - model_options = ''
216   - for key,value in dico[method].iteritems():
217   - model_options = model_options + ' -classifier.{}.{} {}'.format(method, key, value)
218   -
219   - command = command + model_options
220   -
221   - subprocess.call(command, shell=True)
222   -
223   - else:
224   - TrainVectorClassifier = otbApplication.Registry.CreateApplication("TrainVectorClassifier")
225   - TrainVectorClassifier.SetParameterStringList("io.vd", [training_samples_extracted])
226   -
227   - TrainVectorClassifier.SetParameterString("io.stats", img_stats)
228   - TrainVectorClassifier.SetParameterString("io.out", model_out)
229   - TrainVectorClassifier.SetParameterString("classifier", method)
230   -
231   - TrainVectorClassifier.SetParameterStringList("feat", ["band_0", "band_1", "band_2", "band_3"])
232   -
233   - TrainVectorClassifier.UpdateParameters()
234   - TrainVectorClassifier.SetParameterStringList("cfield", [ "class" ])
235   - TrainVectorClassifier.ExecuteAndWriteOutput()
236   -
237   - print('Done')
238   - else:
239   - print("Training not done this time")
240   -
241   - return model_out
242   -
243   -
244   -
245   -
246   -#-------- 3. CLASSIFICATION ---------------------
247   -
248   -def image_classification(main_dir, raw_img, model, img_stats, shell=True, proceed=True, additional_name=''):
249   - '''
250   - 6. Classification on the image
251   - '''
252   - # otbcli_ImageClassifier -in in/cutted.tif -model out/model.rf -out out/labeled_image.tif
253   -
254   - img_labeled = op.join(main_dir, 'Out', "labeled_image{}.tif".format(additional_name))
255   - confidence_map = op.join(main_dir, 'Out', "confidence{}.tif".format(additional_name))
256   -
257   - if proceed==True:
258   - if shell==True:
259   - print(" Image Classification (shell)")
260   - #~ command = 'otbcli_ImageClassifier -in {} -imstat {} -model {} -out {} -confmap {}'.format(raw_img, img_stats, model, img_labeled, confidence_map)
261   - command = 'otbcli_ImageClassifier -in {} -model {} -out {} -confmap {}'.format(raw_img, model, img_labeled, confidence_map)
262   -
263   - subprocess.call(command, shell=True)
264   - else:
265   - print(" Image Classification (API)")
266   - ImageClassifier = otbApplication.Registry.CreateApplication("ImageClassifier")
267   - ImageClassifier.SetParameterString("in", raw_img)
268   - #~ ImageClassifier.SetParameterString("imstat",img_stats)
269   - ImageClassifier.SetParameterString("model", model)
270   - ImageClassifier.SetParameterString("out", img_labeled)
271   - ImageClassifier.SetParameterString("confmap", confidence_map)
272   - ImageClassifier.ExecuteAndWriteOutput()
273   -
274   - print('Done')
275   - else:
276   - print("Classification not done this time")
277   -
278   - return img_labeled
279   -
280   -
281   -def fancy_classif_viz(main_dir, img_labeled, proceed=True):
282   - '''
283   - 7. Fancy visualisation
284   - '''
285   - # otbcli_ColorMapping -in out/labeled_image.tif -method custom -method.custom.lut in/color_table.txt -out out/colorized_classif.png
286   -
287   - color_table = op.join(main_dir, 'Other', 'color_table.txt')
288   - out_image_colorized = op.join(main_dir, 'Out', 'colorized_classif.png')
289   -
290   - ColorMapping = otbApplication.Registry.CreateApplication("ColorMapping")
291   - ColorMapping.SetParameterString("in", img_labeled)
292   - ColorMapping.SetParameterString("method", "custom")
293   - ColorMapping.SetParameterString("method.custom.lut", color_table)
294   - ColorMapping.SetParameterString("out", out_image_colorized)
295   - ColorMapping.UpdateParameters()
296   - ColorMapping.ExecuteAndWriteOutput()
297   -
298   - return out_image_colorized
299   -
300   -
301   -def model_validation(main_dir, img_labeled, testing_shp, show=True, proceed=True):
302   - '''
303   - 8. Confusion matrix
304   - '''
305   - # otbcli_ComputeConfusionMatrix -in out/labeled_image.tif -out out/confusion_matrix.csv -ref vector -ref.vector.in in/testing_clouds.shp -ref.vector.field class
306   -
307   - conf_matrix = op.join(main_dir, 'Statistics', 'confusion_matrix.csv')
308   -
309   - print(' Confusion matrix computing')
310   - ComputeConfusionMatrix = otbApplication.Registry.CreateApplication("ComputeConfusionMatrix")
311   - ComputeConfusionMatrix.SetParameterString("in", img_labeled)
312   - ComputeConfusionMatrix.SetParameterString("ref", "vector")
313   - ComputeConfusionMatrix.SetParameterString("ref.vector.in", testing_shp)
314   - ComputeConfusionMatrix.SetParameterString("out", conf_matrix)
315   - ComputeConfusionMatrix.UpdateParameters()
316   - ComputeConfusionMatrix.SetParameterString("ref.vector.field", "class")
317   - ComputeConfusionMatrix.ExecuteAndWriteOutput()
318   -
319   - print('Done')
320   -
321   - with open(conf_matrix, 'rb') as csvfile:
322   - reader = csv.reader(csvfile, delimiter=',')
323   - rows = []
324   - for row in reader:
325   - rows.append(row)
326   - rows = rows[2:4]
327   - TP = rows[0][0]
328   - FP = rows[0][1]
329   - FN = rows[1][0]
330   - TN = rows[1][1]
331   - print('Confusion matrix:')
332   - print('{:10}'.format(TP) + ' | ' + '{:10}'.format(FP) +
333   - '\n --------------\n' + '{:10}'.format(FN) + ' | ' + '{:10}'.format(TN))
334   -
335   -
336   - return conf_matrix
337   -
338   -
339   -def total_process(main_dir, raw_img, training_shp, testing_shp, method, additional_name=''):
340   - create_directories(main_dir)
341   - training_samples_extracted, img_stats = total_preprocessing(main_dir, raw_img, training_shp)
342   - model = train_model(main_dir, method, training_samples_extracted, img_stats, shell=True, proceed=True)
343   - img_labeled = image_classification(main_dir, raw_img, model, img_stats, shell=True, proceed=True, additional_name = additional_name)
344   - model_validation(main_dir, img_labeled, testing_shp)
345   - fancy_classif_viz(main_dir, img_labeled, proceed=True)
346   -
347   - return img_labeled
348   -
349   -def main():
350   - method = 'rf'
351   -
352   - main_dir = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans'
353   - im_name = 'orleans_bands.tif'
354   - raw_img = op.join(main_dir, 'In_data', 'Image', im_name)
355   - testing_shp = op.join(main_dir, 'In_data', 'Masks', 'validation_points.shp')
356   - training_shp = op.join(main_dir, 'In_data', 'Masks', 'train_points_ext.shp')
357   - #~ training_shp = op.join(main_dir, 'In_data', 'Masks', 'human_reference.shp')
358   -
359   -
360   - create_directories(main_dir)
361   -
362   - #~ img_stats = compute_image_stats(main_dir, raw_img)
363   - #~ print(img_stats)
364   - img_stats = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/Statistics/img_stats.xml'
365   -
366   - proceed = True
367   - class_stats = compute_samples_stats(main_dir, raw_img, training_shp, proceed=True)
368   - training_samples_location = select_samples(main_dir, raw_img, training_shp, class_stats, strategy="constant_8000", proceed=proceed)
369   - training_samples_extracted = extract_samples(main_dir, raw_img, training_samples_location, proceed=proceed)
370   -
371   - model = train_model(main_dir, method, training_samples_extracted, img_stats, shell=True, proceed=True)
372   - additional_name = '0'
373   - img_labeled = image_classification(main_dir, raw_img, model, img_stats, shell=True, proceed=True, additional_name = additional_name)
374   - model_validation(main_dir, img_labeled, testing_shp)
375   - fancy_classif_viz(main_dir, img_labeled, proceed=True)
376   -
377   -
378   - #~ # Step 1 : human labellisation
379   - #~ training_shp = op.join(main_dir, 'In_data', 'Masks', 'training_clouds2.shp')
380   - #~ img_labeled = total_process(main_dir, raw_img, training_shp, testing_shp, method, additional_name='1')
381   -
382   - #~ img_labeled = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/Out/labeled_image1.tif'
383   -
384   - #~ # Step 2 : ML labellisation
385   - #~ training_shp = op.join(main_dir, 'Intermediate', 'RF_clouds1.shp')
386   - #~ resample(img_labeled, training_shp)
387   - #~ img_labeled = total_process(main_dir, raw_img, training_shp, testing_shp, method, additional_name='2')
388   -
389   - #~ # Step 3 : ML labellisation
390   - #~ training_shp = op.join(main_dir, 'Intermediate', 'RF_clouds2.shp')
391   - #~ resample(img_labeled, training_shp)
392   - #~ img_labeled = total_process(main_dir, raw_img, training_shp, testing_shp, method, additional_name='3')
393   -
394   - #~ # Step 4 : ML labellisation
395   - #~ training_shp = op.join(main_dir, 'Intermediate', 'RF_clouds3.shp')
396   - #~ resample(img_labeled, training_shp)
397   - #~ img_labeled = total_process(main_dir, raw_img, training_shp, testing_shp, method, additional_name='4')
398   -
399   -
400   - #~ training_shp = op.join(main_dir, 'Intermediate', 'RF_clouds1.shp')
401   - #~ resample(img_labeled, training_shp)
402   -
403   - #~ img_stats = compute_image_stats(main_dir, raw_img)
404   -
405   -
406   - #~ training_samples_extracted = samples_preprocessing(main_dir, raw_img, training_shp, proceed=True)
407   -
408   - #~ model = train_model(main_dir, method, training_samples_extracted, img_stats, shell=True, proceed=True)
409   - #~ img_labeled = image_classification(main_dir, raw_img, model, img_stats, shell=True, proceed=True, additional_name = '2')
410   - #~ model_validation(main_dir, img_labeled, testing_shp)
411   - #~ fancy_classif_viz(main_dir, img_labeled, proceed=True)
412   -
413   -
414   -
415   -if __name__ == '__main__':
416   - main()
417   -
418   -
419   -
420   -
421   -
ALCD/unused/PCA_analysis.py
... ... @@ -1,95 +0,0 @@
1   -#!/usr/bin/python
2   -
3   -
4   -#~ import numpy as np
5   -#~ from matplotlib.mlab import PCA
6   -
7   -#~ data = np.array(np.random.randint(10,size=(10,3)))
8   -#~ results = PCA(data)
9   -
10   -#~ print(results)
11   -import sqlite3
12   -import matplotlib.pyplot as plt
13   -import numpy as np
14   -from matplotlib.mlab import PCA
15   -import otbApplication
16   -import os.path as op
17   -
18   -
19   -def parse_sqlite(sqlite_file):
20   - '''Parse a sqlite file, and return the headers along with the data
21   - '''
22   - conn = sqlite3.connect(sqlite_file)
23   - cur = conn.cursor()
24   - data = cur.execute("SELECT * FROM output")
25   -
26   - # Get the headers of the columns
27   - headers = data.description
28   - headers = [h[0] for h in headers]
29   -
30   - # Get all the data
31   - data_all = data.fetchall()
32   -
33   - return headers, data_all
34   -
35   -def plot_histogram(data):
36   - N_points = 100000
37   - n_bins = 200
38   -
39   - fig, axs = plt.subplots(1, 2, sharey=True, tight_layout=True)
40   -
41   - # We can set the number of bins with the `bins` kwarg
42   - axs[0].hist(data, bins=n_bins)
43   - axs[1].hist(data, bins=n_bins)
44   -
45   - plt.show()
46   -
47   -def PCA_analysis(data):
48   - results = PCA(data)
49   -
50   - print(results)
51   -
52   -def dimensionality_reduction(in_tif, out_tif):
53   - #!/usr/bin/python
54   -
55   - # Import the otb applications package
56   - import otbApplication
57   -
58   - # The following line creates an instance of the DimensionalityReduction application
59   - DimensionalityReduction = otbApplication.Registry.CreateApplication("DimensionalityReduction")
60   -
61   - # The following lines set all the application parameters:
62   - DimensionalityReduction.SetParameterString("in", in_tif)
63   -
64   - DimensionalityReduction.SetParameterString("out", out_tif)
65   - DimensionalityReduction.SetParameterString("nbcomp", "1")
66   - DimensionalityReduction.SetParameterString("normalize", "true")
67   -
68   - DimensionalityReduction.SetParameterString("method","pca")
69   - DimensionalityReduction.UpdateParameters()
70   -
71   - # The following line execute the application
72   - DimensionalityReduction.ExecuteAndWriteOutput()
73   -
74   -
75   -def main():
76   - #~ sqlite_file = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/Samples/training_samples_extracted.sqlite'
77   - #~ headers, data_all = parse_sqlite(sqlite_file)
78   - #~ print(headers)
79   - #~ print(len(data_all))
80   -
81   - #~ band_0_idx = headers.index('band_0')
82   - #~ print(band_0_idx)
83   - #~ band0 = [d[band_0_idx] for d in data_all]
84   -
85   - #~ plot_histogram(band0)
86   -
87   - mdir = '/mnt/data/home/baetensl/OTB_codes/OTB_commands/Full_orleans/In_data/Image'
88   - in_tif = op.join(mdir, "all_bands.tif")
89   - out_tif = op.join(mdir, "reduced.tif")
90   - dimensionality_reduction(in_tif, out_tif)
91   -
92   -
93   -
94   -if __name__ == '__main__':
95   - main()
ALCD/unused/QGIS_applystyle.py
... ... @@ -1,28 +0,0 @@
1   -# -*- coding: utf-8 -*-
2   -"""
3   -Created on Mon Apr 23 13:39:10 2018
4   -
5   -@author: baetensl
6   -"""
7   -
8   -from qgis.core import *
9   -import os
10   -#copy line 9-21 and change file names and group names if you have more groups
11   -
12   -QML_file = ('/mnt/data/home/baetensl/OTB_codes/style_RVB.qml')#insert path to qml file
13   -#add other qml files if you want to change style for more groups
14   -
15   -
16   -def applystyle_group(name):
17   - root = QgsProject.instance().layerTreeRoot()
18   - point = root.findGroup(name) #Find Group
19   - for child in point.children():
20   - if isinstance(child, QgsLayerTreeLayer):
21   - if child.layer().type()==0:
22   - child.layer().loadNamedStyle(QML_file)#change the file name accordingly
23   - #you can add styles for other types of layers in the same group (line, point and polygon)
24   -
25   -try: #If group is not present this will keep script running if you want to add more
26   - applystyle_group("raster")#insert name of QGIS group
27   -except Exception:
28   - pass
29 0 \ No newline at end of file
ALCD/unused/QGIS_launch.py
... ... @@ -1,79 +0,0 @@
1   -# -*- coding: utf-8 -*-
2   -"""
3   -Created on Mon Apr 23 11:36:13 2018
4   -
5   -@author: baetensl
6   -"""
7   -
8   -import os
9   -
10   -def compose_bands(red_band_num=4, green_band_num=3, blue_band_num=2, additional_band_num = 10):
11   - '''
12   - Compose the given bands, and return 3 bands only, so that the user can
13   - visualize a specific set of composition
14   - '''
15   -
16   - L1C_dir = '/mnt/data/SENTINEL2/L1C_PDGS/'
17   - sub_dir = 'Orleans/S2A_MSIL1C_20170615T105031_N0205_R051_T31UDP_20170615T105505.SAFE/GRANULE/L1C_T31UDP_A010344_20170615T105505/IMG_DATA/'
18   -
19   - tile = 'T31UDP'
20   - date = '20170615'
21   -
22   - common_file = 'T31UDP_20170615T105031_B'#05.jp2
23   -
24   - in_dir = L1C_dir + sub_dir
25   -
26   - print(in_dir)
27   -
28   - # to create the composition of multiple bands into one tiff
29   - current_dir = '/mnt/data/home/baetensl/OTB_codes'
30   - out_dir = os.path.join(current_dir, 'OTB_in') #TODO : faire la creation du folder automatiquement
31   -
32   - out_name = 'B{:02d}_B{:02d}_B{:02d}_B{:02d}.tif'.format(red_band_num, green_band_num, blue_band_num, additional_band_num)
33   - outTIF = os.path.join(out_dir, out_name) #output tiff
34   -
35   - selected_bands = [red_band_num, green_band_num, blue_band_num, additional_band_num]
36   - red_band = os.path.join(in_dir, (common_file+'{num:02d}'.format(num=selected_bands[0])+'.jp2'))
37   - green_band = os.path.join(in_dir, (common_file+'{num:02d}'.format(num=selected_bands[1])+'.jp2'))
38   - blue_band = os.path.join(in_dir, (common_file+'{num:02d}'.format(num=selected_bands[2])+'.jp2'))
39   - additional_band = os.path.join(in_dir, (common_file+'{num:02d}'.format(num=selected_bands[3])+'.jp2'))
40   -
41   - build_vrt = 'gdalbuildvrt -separate {} {} {} {} {}'.format(outTIF, red_band, green_band, blue_band, additional_band)
42   - print(build_vrt)
43   - os.system(build_vrt)
44   -
45   - return outTIF
46   -
47   -def open_in_QGIS(filename):
48   - pythonCode = '/mnt/data/home/baetensl/OTB_codes/QGIS_applystyle.py'
49   - launch_qgis = 'qgis {} --code {}'.format(filename, pythonCode)
50   - launch_qgis = 'qgis {}'.format(filename)
51   -
52   - print(launch_qgis)
53   -# launch_qgis = 'qgis --noplugins {}'.format(filename)
54   -
55   - os.system(launch_qgis)
56   -
57   -
58   -
59   -#main_dir = '/mnt/data/home/baetensl/GDAL/'
60   -#source_raster = main_dir + "MAJA_data/S2A_OPER_SSC_PDTIMG_L2VALD_31UDP____20161217_FRE_R1.DBL.TIF"
61   -#
62   -## save all the bands in separate files
63   -#for band in [1,2,3]:
64   -# destination_raster = main_dir + "python_out/band{}.tif".format(band)
65   -# command = 'gdal_translate "{}" -b {} "{}"'.format(source_raster, band, destination_raster)
66   -# print(command)
67   -## os.system(command)
68   -
69   -
70   -def main():
71   - outFile = compose_bands(4,3,2,10)
72   - open_in_QGIS(outFile)
73   -
74   -if __name__ == '__main__':
75   - main()
76   -
77   -
78   -
79   -
ALCD/unused/features_visualization.py
... ... @@ -1,331 +0,0 @@
1   -#!/usr/bin/python
2   -# -*- coding: utf-8 -*-
3   -
4   -import sys
5   -import os
6   -import os.path as op
7   -import glob
8   -import random
9   -import numpy as np
10   -import matplotlib.pyplot as plt
11   -import matplotlib
12   -import json
13   -import otbApplication
14   -import csv
15   -
16   -# TODO 13/07/2018 : ADD COMMENTS IN THIS CODE
17   -
18   -def extract_class_band(raw_img, img_labeled, class_nb, band):
19   - '''
20   - From the in_tif and the labeled image, extract the band for the given
21   - class
22   - Mostly used to plot histograms for each class and band
23   - '''
24   -
25   - BandMathX = otbApplication.Registry.CreateApplication("BandMathX")
26   - BandMathX.SetParameterStringList("il", [str(raw_img), str(img_labeled)])
27   -
28   - if str(class_nb) == "-1":
29   - # No class is negative, so it will take the full image
30   - expression = "(im2b1 != -1 ? im1b{band}+0 : 0)".format(band = str(band))
31   - else:
32   - # For all the other classes, will extract the appropriate pixels
33   - expression = "(im2b1 == {class_nb} ? im1b{band} : 0)".format(class_nb = str(class_nb), band = str(band))
34   -
35   - BandMathX.SetParameterString("exp", expression)
36   - BandMathX.UpdateParameters()
37   - BandMathX.Execute()
38   -
39   - # Pass directly to numpy array, doesnt write on disk
40   - extraction_output = BandMathX.GetImageAsNumpyArray('out')
41   - data = np.array(extraction_output)
42   -
43   - return data
44   -
45   -
46   -
47   -
48   -def compute_correlation(global_parameters, class_nb, features):
49   - '''
50   - Compute the correlation between N features for a given class
51   - '''
52   - main_dir = global_parameters["user_choices"]["main_dir"]
53   - raw_img = op.join(main_dir, 'In_data', 'Image', global_parameters["user_choices"]["raw_img"])
54   - img_labeled = op.join(main_dir, 'Out', global_parameters["general"]["img_labeled_regularized"])
55   -
56   - # Load the pixel values for each feature of the image
57   - vectors = []
58   - valid_index = []
59   - invalid_indexes = []
60   - for feat in features:
61   - data = extract_class_band(raw_img, img_labeled, class_nb, feat)
62   - vector_temp = np.ravel(data)
63   - # Remove the NaN values if they exist (e.g. for the NDVI and NDWI)
64   - vector_temp[np.isnan(vector_temp)] = 0
65   -
66   - # A pixel is valid if its not 0 (i.e. doesnt belong to the class)
67   - # it is needed to remove all 0 values (otherwise it increases the correlation
68   - # artificially)
69   - if valid_index == []:
70   - valid_index = [vector_temp != 0]
71   - else:
72   - valid_index = valid_index and [vector_temp != 0]
73   -
74   - if len(vector_temp[valid_index]) == 0:
75   - print('No valid data for the class {}, correlation not computed'.format(class_nb))
76   - return
77   -
78   - vectors.append(vector_temp)
79   - # Takes only the valid indices
80   - vectors = [vec[valid_index] for vec in vectors]
81   -
82   - # Compute the correlation between features
83   - correlation = np.corrcoef(vectors)
84   - correlation_float = np.copy(correlation)
85   -
86   - # Format to export nicely to a .csv
87   - for row in correlation:
88   - for j in range(len(row)):
89   - row[j] = '{0:.2f}'.format(row[j])
90   -
91   - features = ['feat {}'.format(f) for f in features]
92   - correlation = np.row_stack([features, correlation])
93   -
94   - features.insert(0,'class {}'.format(class_nb))
95   - correlation = np.column_stack((features, correlation))
96   -
97   - correlation_csv = op.join(main_dir, 'Statistics',
98   - 'correlations', 'corr_table_{}.csv'.format(class_nb))
99   -
100   - np.savetxt(correlation_csv, correlation, delimiter=',', fmt = '%s')
101   - features = features[1:]
102   -
103   - # Display a matrix of correlation with matplotlib
104   - display_correlation(global_parameters, features, correlation_float, class_nb)
105   -
106   - return features, correlation_float
107   -
108   -
109   -
110   -def display_correlation(global_parameters, features, correlation, class_nb):
111   - '''
112   - Plot the correlation matrix between different features
113   - '''
114   - # Get the class names and numbers
115   - all_class_names = []
116   - all_class_nb = []
117   - all_class_names.append("full_image")
118   - all_class_nb.append(-1)
119   - masks = global_parameters["masks"]
120   - for m in masks:
121   - all_class_names.append(str(m))
122   - all_class_nb.append(int(global_parameters["masks"][m]["class"]))
123   -
124   - try:
125   - class_name = all_class_names[all_class_nb.index(class_nb)]
126   - except:
127   - class_name = 'null_value'
128   -
129   - # Rename the features for better readability
130   - features = [('f' + f.split(' ')[1]) for f in features]
131   -
132   - # Takes the absolute value as 1 and -1 are well correlated
133   - correlation = np.absolute(correlation)
134   -
135   - fig = plt.figure()
136   - ax = fig.add_subplot(111)
137   - cax = ax.matshow(correlation, interpolation='nearest', vmin=0, vmax=1, cmap='winter')
138   - fig.colorbar(cax)
139   -
140   - # To have the correct labels
141   - ax.set_xticklabels(features)
142   - ax.set_yticklabels(features)
143   - ax.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(1))
144   - ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(1))
145   - plt.title('abs(correlation) for {}'.format(class_name))
146   - plt.show(block = False)
147   -
148   - # Save the figure