Commit d49313149a59dfbcd7a8fc86d5a8f020b6e61441

Authored by Louis Baetens
1 parent 1de4b8b9
Exists in master

ENH pixels features op command line possible

PCC/parameters_files/comparison_parameters.json
... ... @@ -7,13 +7,13 @@
7 7 "resolution": "60"
8 8 },
9 9 "automatically_generated": {
10   - "alcd_main_dir": "/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/Ispra_32TMR_20170815",
11   - "current_date": "20170815",
12   - "location": "Ispra",
13   - "main_dir": "/mnt/data/home/baetensl/clouds_detection_git/Data_PCC/HOT016/Ispra_32TMR_20170815",
14   - "raw_img": "Ispra_bands.tif",
  10 + "alcd_main_dir": "/mnt/data/home/baetensl/clouds_detection_git/Data_ALCD/Arles_31TFJ_20171002",
  11 + "current_date": "20171002",
  12 + "location": "Arles",
  13 + "main_dir": "/mnt/data/home/baetensl/clouds_detection_git/Data_PCC/HOT016/Arles_31TFJ_20171002",
  14 + "raw_img": "Arles_bands.tif",
15 15 "sub_data_dir": "HOT016",
16   - "tile": "32TMR"
  16 + "tile": "31TFJ"
17 17 },
18 18 "maja_parameters": {
19 19 "maja_dir_name": "MAJA_3_1_S2AS2B_HOT016",
... ...
PCC/pixels_features_analysis.py
... ... @@ -9,87 +9,16 @@ import json
9 9 import otbApplication
10 10 from matplotlib.lines import Line2D
11 11 import collections
  12 +import argparse
12 13  
13   -#~ def extract_class_band(raw_img, band_to_extract, diff_img, diff_class_nb, alcd_labeled_img='', alcd_class_nb = 0):
14   - #~ '''
15   - #~ From the difference image (e.g. diff_binary) between the processing
16   - #~ chain and the ALCD, and the ALCD labeled, extract the data for the pixels
17   - #~ in the class_nb for both images
18   - #~ Mostly used to plot histograms for each class and band
19   - #~ '''
20   -
21   - #~ BandMathX = otbApplication.Registry.CreateApplication("BandMathX")
22   - #~ if alcd_labeled_img == '':
23   - #~ BandMathX.SetParameterStringList("il", [str(raw_img), str(diff_img)])
24   - #~ expression = "(im2b1 == {diff_class_nb} ? im1b{band} : 0)".format(
25   - #~ diff_class_nb = str(diff_class_nb), band = str(band_to_extract))
26   - #~ else:
27   - #~ BandMathX.SetParameterStringList("il", [str(raw_img), str(diff_img), str(alcd_labeled_img)])
28   - #~ expression = "((im2b1 == {diff_class_nb} && im3b1 == {alcd_class_nb}) ? im1b{band} : 0)".format(
29   - #~ diff_class_nb = str(diff_class_nb),
30   - #~ alcd_class_nb = str(alcd_class_nb),
31   - #~ band = str(band_to_extract))
32   -
33   - #~ BandMathX.SetParameterString("exp", expression)
34   -
35   - #~ BandMathX.UpdateParameters()
36   - #~ BandMathX.Execute()
37   -
38   - #~ # Pass directly to numpy array, doesnt write on disk
39   - #~ extraction_output = BandMathX.GetImageAsNumpyArray('out')
40   - #~ data = np.array(extraction_output)
41   -
42   - #~ return data
43   -
44   -
45   -
46   -#~ def count_intersected_classes(diff_img, diff_class_nb, alcd_labeled_img, alcd_class_nb):
47   - #~ '''
48   - #~ Count the number of pixels being at the same time of the
49   - #~ diff_class_nb and alcd_class_nb.
50   - #~ e.g. if diff_class_nb = 1 and alcd_class_nb = 2, it will returns the
51   - #~ number of true positive and low_clouds pixels.
52   - #~ '''
53   -
54   - #~ Superimpose = otbApplication.Registry.CreateApplication("Superimpose")
55   - #~ Superimpose.SetParameterString("inr", str(diff_img))
56   - #~ Superimpose.SetParameterString("inm", str(alcd_labeled_img))
57   - #~ Superimpose.SetParameterString("interpolator", "nn")
58   - #~ Superimpose.UpdateParameters()
59   - #~ Superimpose.Execute()
60   -
61   -
62   - #~ BandMathX = otbApplication.Registry.CreateApplication("BandMathX")
63   -
64   - #~ BandMathX.AddParameterStringList("il", str(diff_img))
65   - #~ BandMathX.AddImageToParameterInputImageList("il", Superimpose.GetParameterOutputImage("out"))
66   - #~ #im1 is the diff image, im2 the ALCD ref
67   - #~ expression = "((im1b1 == {diff_class_nb} && im2b1 == {alcd_class_nb}) ? im1b1 : 0)".format(
68   - #~ diff_class_nb = str(diff_class_nb),
69   - #~ alcd_class_nb = str(alcd_class_nb))
70   -
71   - #~ BandMathX.SetParameterString("exp", expression)
72   -
73   - #~ BandMathX.UpdateParameters()
74   - #~ BandMathX.Execute()
75   -
76   - #~ # Pass directly to numpy array, doesnt write on disk
77   - #~ extraction_output = BandMathX.GetImageAsNumpyArray('out')
78   - #~ data = np.array(extraction_output)
79   -
80   - #~ return data
81   -
82   -
83   -
84   -
85   -
86   -
  14 +
87 15  
88 16 def get_alcd_label_from_diff(diff_img, diff_class_nb, alcd_labeled_img):
89 17 '''
90   - A FAIRE
  18 + Get the pixels of the ALCD labeled image, which are of the class
  19 + diff_class_nb (linked to the confusion matrix, so 1 2 3 or 4) of
  20 + the image diff_img
91 21 '''
92   -
93 22 Superimpose = otbApplication.Registry.CreateApplication("Superimpose")
94 23 Superimpose.SetParameterString("inr", str(diff_img))
95 24 Superimpose.SetParameterString("inm", str(alcd_labeled_img))
... ... @@ -103,7 +32,9 @@ def get_alcd_label_from_diff(diff_img, diff_class_nb, alcd_labeled_img):
103 32 BandMathX.AddParameterStringList("il", str(diff_img))
104 33 BandMathX.AddImageToParameterInputImageList("il", Superimpose.GetParameterOutputImage("out"))
105 34 #im1 is the diff image, im2 the ALCD ref
106   - #~ expression = "(im1b1 == {diff_class_nb} ? im2b1 : -1)".format(
  35 + # pixels that have a null value / no-data will be set to 0,
  36 + # the tones dont belong to the class will be set to -1,
  37 + # and the others to their labeled class
107 38 expression = "((im1b1 == {diff_class_nb} ? im2b1 : -1) + (im1b1 == 0 ? 1 : 0))".format(
108 39 diff_class_nb = str(diff_class_nb))
109 40  
... ... @@ -116,8 +47,10 @@ def get_alcd_label_from_diff(diff_img, diff_class_nb, alcd_labeled_img):
116 47 extraction_output = BandMathX.GetImageAsNumpyArray('out')
117 48 alcd_masked = np.array(extraction_output).flatten().astype(int)
118 49  
  50 + # dict with number of pixels per class
119 51 class_numbers = collections.Counter(alcd_masked)
120   -
  52 +
  53 + # number of valid pixels on the image
121 54 nb_valid_pixels = len(alcd_masked)-class_numbers[0]
122 55  
123 56  
... ... @@ -126,8 +59,75 @@ def get_alcd_label_from_diff(diff_img, diff_class_nb, alcd_labeled_img):
126 59  
127 60  
128 61  
  62 +def get_miclassified_histograms(location, date, cas_alcd, processing_chain, pcc_sub_dir = None, diff_class_nb = 1, display = False):
  63 + paths_configuration = json.load(open(op.join('..', 'paths_configuration.json')))
  64 +
  65 + data_alcd_path = paths_configuration["data_paths"]["data_alcd"]
  66 + data_pcc_path = paths_configuration["data_paths"]["data_pcc"]
  67 + if pcc_sub_dir != None:
  68 + data_pcc_path = op.join(data_pcc_path, pcc_sub_dir)
  69 + tile = paths_configuration["tile_location"][location]
  70 +
  71 + location_dir = location + '_' + tile + '_' + date
  72 + raw_img = glob.glob(op.join(data_alcd_path, location_dir, 'In_data', 'Image', (location + '_bands.tif')))[0]
  73 +
  74 + main_dir = op.join(data_pcc_path, '{}_{}_{}'.format(location, tile, date))
  75 +
  76 + alcd_labeled_img = op.join(main_dir, 'Original_data', 'labeled_img_original.tif')
  77 +
  78 + if cas_alcd == 1:
  79 + reference_sub_dir = 'ALCD_initial'
  80 + if processing_chain == 'maja':
  81 + processing_name = 'maja_ini'
  82 + elif processing_chain == 'sen2cor':
  83 + processing_name = 'sen2cor_ini'
  84 + elif processing_chain == 'fmask':
  85 + processing_name = 'fmask_ini'
  86 + elif cas_alcd == 2:
  87 + reference_sub_dir = 'ALCD_dilated'
  88 + if processing_chain == 'maja':
  89 + processing_name = 'maja_ini'
  90 + elif processing_chain == 'sen2cor':
  91 + processing_name = 'sen2cor_dilat'
  92 + elif processing_chain == 'fmask':
  93 + processing_name = 'fmask_dilat'
  94 + elif cas_alcd == 3:
  95 + reference_sub_dir = 'ALCD_initial'
  96 + if processing_chain == 'maja':
  97 + processing_name = 'maja_erode'
  98 + elif processing_chain == 'sen2cor':
  99 + processing_name = 'sen2cor_ini'
  100 + elif processing_chain == 'fmask':
  101 + processing_name = 'fmask_ini'
  102 +
  103 +
  104 +
  105 + diff_img = glob.glob(op.join(main_dir, 'Binary_difference', reference_sub_dir, '*{}*.tif'.format(processing_name)))[0]
  106 +
  107 + numbers_dict, nb_valid_pixels = get_alcd_label_from_diff(diff_img, diff_class_nb, alcd_labeled_img)
  108 +
  109 +
  110 + non_zero_class_nb = []
  111 + nb_per_class = []
  112 +
  113 +
  114 + all_classes_nb = np.arange(1,8)
  115 + for alcd_class_nb in all_classes_nb:
  116 + if numbers_dict[alcd_class_nb] > 0:
  117 + non_zero_class_nb.append(alcd_class_nb)
  118 + nb_per_class.append(numbers_dict[alcd_class_nb])
  119 +
  120 +
  121 + total_samples = sum(nb_per_class)
  122 + nb_per_class = [float(n)/total_samples for n in nb_per_class]
  123 +
  124 + proportion_in_image = float(total_samples)/float(nb_valid_pixels)
  125 +
  126 + return non_zero_class_nb, nb_per_class, proportion_in_image
129 127  
130   -def plot_confusion_repartition(location, date, processing_chain, cas_alcd):
  128 +
  129 +
  130 +def plot_confusion_repartition(location, date, processing_chain, cas_alcd, pcc_sub_dir = None):
131 131 '''
132 132 Plot the repartition of classes for each cell of the confusion
133 133 matrix.
... ... @@ -135,6 +135,9 @@ def plot_confusion_repartition(location, date, processing_chain, cas_alcd):
135 135 '''
136 136 paths_configuration = json.load(open(op.join('..', 'paths_configuration.json')))
137 137 data_pcc_dir = paths_configuration["data_paths"]["data_pcc"]
  138 + if pcc_sub_dir != None:
  139 + data_pcc_dir = op.join(data_pcc_dir, pcc_sub_dir)
  140 +
138 141 tile = paths_configuration["tile_location"][location]
139 142 main_dir = op.join(data_pcc_dir, '{}_{}_{}'.format(location, tile, date))
140 143  
... ... @@ -155,7 +158,8 @@ def plot_confusion_repartition(location, date, processing_chain, cas_alcd):
155 158  
156 159 for confusion_mat_place in [1,2,3,4]:
157 160 non_zero_class_nb, nb_per_class, proportion_in_image = get_miclassified_histograms(location,
158   - date, cas_alcd, processing_chain, diff_class_nb = confusion_mat_place, display = True)
  161 + date, cas_alcd, processing_chain, pcc_sub_dir = pcc_sub_dir,
  162 + diff_class_nb = confusion_mat_place, display = True)
159 163 class_nb.append(non_zero_class_nb)
160 164 all_proportion.append(nb_per_class)
161 165 mat_conf_proportion.append(proportion_in_image)
... ... @@ -209,131 +213,11 @@ def plot_confusion_repartition(location, date, processing_chain, cas_alcd):
209 213 plt.tight_layout()
210 214  
211 215 # save the figure
212   - plt.show()
213   - out_fig = op.join(main_dir, 'Statistics', 'features', 'repartition_matconf_{}_{}_{}.png'.format(processing_chain, location, date))
214   - #~ plt.savefig(out_fig, bbox_inches='tight')
215   -
216   -
217   -
218   -
219   -
220   -
221   -
222   -
223   -
  216 + #~ plt.show()
  217 + out_fig = op.join(main_dir, 'Statistics', 'features', 'repartition_matconf_{}_{}_{}_case{}.png'.format(processing_chain, location, date, cas_alcd))
  218 + plt.savefig(out_fig, bbox_inches='tight')
224 219  
225   -
226   -def get_miclassified_histograms(location, date, cas_alcd, processing_chain, diff_class_nb = 1, display = False):
227   - paths_configuration = json.load(open(op.join('..', 'paths_configuration.json')))
228   -
229   - data_alcd_path = paths_configuration["data_paths"]["data_alcd"]
230   - data_pcc_path = paths_configuration["data_paths"]["data_pcc"]
231   - tile = paths_configuration["tile_location"][location]
232   -
233   - location_dir = location + '_' + tile + '_' + date
234   - raw_img = glob.glob(op.join(data_alcd_path, location_dir, 'In_data', 'Image', (location + '_bands.tif')))[0]
235   -
236   - main_dir = op.join(data_pcc_path, '{}_{}_{}'.format(location, tile, date))
237   -
238   - alcd_labeled_img = op.join(main_dir, 'Original_data', 'labeled_img_original.tif')
239   -
240   - if cas_alcd == 1:
241   - reference_sub_dir = 'ALCD_initial'
242   - if processing_chain == 'maja':
243   - processing_name = 'maja_ini'
244   - elif processing_chain == 'sen2cor':
245   - processing_name = 'sen2cor_ini'
246   - elif processing_chain == 'fmask':
247   - processing_name = 'fmask_ini'
248   - elif cas_alcd == 2:
249   - reference_sub_dir = 'ALCD_dilat'
250   - if processing_chain == 'maja':
251   - processing_name = 'maja_ini'
252   - elif processing_chain == 'sen2cor':
253   - processing_name = 'sen2cor_dilat'
254   - elif processing_chain == 'fmask':
255   - processing_name = 'fmask_dilat'
256   - elif cas_alcd == 3:
257   - reference_sub_dir = 'ALCD_initial'
258   - if processing_chain == 'maja':
259   - processing_name = 'maja_erode'
260   - elif processing_chain == 'sen2cor':
261   - processing_name = 'sen2cor_ini'
262   - elif processing_chain == 'fmask':
263   - processing_name = 'fmask_ini'
264   -
265   -
266   -
267   -
268   - diff_img = glob.glob(op.join(main_dir, 'Binary_difference', reference_sub_dir, '*{}*.tif'.format(processing_name)))[0]
269   -
270   - numbers_dict, nb_valid_pixels = get_alcd_label_from_diff(diff_img, diff_class_nb, alcd_labeled_img)
271   -
272   -
273   - non_zero_class_nb = []
274   - nb_per_class = []
275   -
276   -
277   - all_classes_nb = np.arange(1,8)
278   - for alcd_class_nb in all_classes_nb:
279   - if numbers_dict[alcd_class_nb] > 0:
280   - non_zero_class_nb.append(alcd_class_nb)
281   - nb_per_class.append(numbers_dict[alcd_class_nb])
282   -
283   -
284   - total_samples = sum(nb_per_class)
285   - nb_per_class = [float(n)/total_samples for n in nb_per_class]
286 220  
287   - proportion_in_image = float(total_samples)/float(nb_valid_pixels)
288   -
289   - return non_zero_class_nb, nb_per_class, proportion_in_image
290   -
291   -
292   -
293   -
294   - alcd_class_dict = {'0':'null value', '1':'background', '2':'low clouds',
295   - '3':'high clouds', '4':'clouds shadows', '5':'land',
296   - '6':'water', '7':'snow'}
297   -
298   - # Create the colors
299   - all_classes_nb = np.arange(1,8)
300   - colormap = create_colormap(color_pattern = 'labeled', class_call = 'number')
301   -
302   - non_zero_names = []
303   - non_zero_colors = []
304   - non_zero_class_nb = []
305   -
306   - # get the data for each class
307   - nb_per_class = []
308   - for alcd_class_nb in all_classes_nb:
309   - #~ data = count_intersected_classes(diff_img, diff_class_nb, alcd_labeled_img, alcd_class_nb)
310   - data = get_alcd_label_from_diff(diff_img, diff_class_nb, alcd_labeled_img)
311   - return
312   - data = data.flatten()
313   - total_pixels_in_img = len(data)
314   - # remove the zero values
315   - data = data[data != 0]
316   - data = data[np.isnan(data) == False]
317   -
318   - nb_of_samples = len(data)
319   - if nb_of_samples != 0:
320   - nb_per_class.append(nb_of_samples)
321   - non_zero_names.append(alcd_class_dict[str(alcd_class_nb)])
322   - non_zero_colors.append(colormap[str(alcd_class_nb)])
323   - non_zero_class_nb.append(alcd_class_nb)
324   -
325   -
326   - total_samples = sum(nb_per_class)
327   - nb_per_class = [float(n)/total_samples for n in nb_per_class]
328   -
329   - proportion_in_image = float(total_samples)/float(total_pixels_in_img)
330   -
331   -
332   -
333   - return non_zero_class_nb, nb_per_class, proportion_in_image
334   -
335   -
336   -
337 221  
338 222 def create_colormap(color_pattern = 'diff', class_call = 'number'):
339 223 '''
... ... @@ -383,26 +267,30 @@ def create_colormap(color_pattern = 'diff', class_call = 'number'):
383 267  
384 268  
385 269 def main():
386   - #~ comparison_parameters = json.load(open(op.join('parameters_files','comparison_parameters.json')))
387   -
388   - location = 'Arles'
389   - date = '20171002'
390   - processing_chain = 'maja'
391   - cas_alcd = 1
392   - plot_confusion_repartition(location, date, processing_chain, cas_alcd)
393   -
394   - return
395   -
396   -
397   - for processing_chain in ['maja', 'sen2cor', 'fmask']:
398   - try:
399   - plot_confusion_repartition(comparison_parameters, processing_chain)
400   - except:
401   - pass
402   -
403   - return
404   - alcd_class_nb = 2
405   - extract_bands_per_class(comparison_parameters, processing_chain, alcd_class_nb)
  270 + parser = argparse.ArgumentParser()
  271 +
  272 + parser.add_argument('-l', action='store', default=None, dest='location', help='Name of the location', required=True)
  273 + parser.add_argument('-d', action='store', default=None, dest='date', help='Date to perform', required=True)
  274 + parser.add_argument('-s', action='store', default=None, dest='pcc_sub_dir', help='Name of the PCC subdir (e.g. HOT016)')
  275 + parser.add_argument('-chain', action='store', default=None, dest='processing_chain', help='Name of the processing_chain (maja, sen2cor or fmask)')
  276 + parser.add_argument('-case', action='store', default=2, dest='cas_alcd', help='Number of the ALCD case (1, 2 or 3)')
  277 + results = parser.parse_args()
  278 +
  279 + pcc_sub_dir = results.pcc_sub_dir
  280 + cas_alcd = int(results.cas_alcd)
  281 + location = results.location
  282 + date = results.date
  283 + processing_chain = results.processing_chain
  284 +
  285 + #~ location = 'Arles'
  286 + #~ date = '20171002'
  287 + #~ processing_chain = 'fmask'
  288 + #~ cas_alcd = 1
  289 +
  290 + plot_confusion_repartition(location, date, processing_chain, cas_alcd, pcc_sub_dir)
  291 +
  292 +
  293 +
406 294  
407 295  
408 296  
... ...