Skip to content

s1_cslc_qa

Class to compute stats for geocoded raster and corrections

QualityAssuranceCSLC

Class to compute stats for geocoded raster and corrections

Source code in src/compass/s1_cslc_qa.py 
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
class QualityAssuranceCSLC:
    '''
    Class to compute stats for geocoded raster and corrections
    '''
    stat_names = ['mean', 'min', 'max', 'std']

    def __init__(self):
        self.stats_dict = {}
        self.pixel_percentage_dict = {}
        self.rfi_dict = {}
        self.is_safe_corrupt = False
        self.orbit_dict = {}
        self.output_to_json = False


    def compute_CSLC_raster_stats(self, cslc_h5py_root, bursts):
        '''
        Compute CSLC raster stats. Stats written to HDF5 and saved to class
        dict for later JSON output

        Parameters
        ----------
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5
        bursts: list
            Bursts whose geocoded raster stats are to be computed
        '''
        for b in bursts:
            pol = b.polarization

            # get dataset and compute stats according to dtype
            pol_path = f'{DATA_PATH}/{pol}'
            pol_arr = cslc_h5py_root[pol_path][()]

            # create dict for current polarization
            self.stats_dict[pol] = {}
            pol_dict = self.stats_dict[pol]

            # compute power or phase then write stats to HDF5 for CSLC
            for pwr_phase in ['power', 'phase']:
                # create dict to store real/imaginary stat items
                pol_dict[pwr_phase] = {}

                # create HDF5 group for power or phase stats of current
                # polarization
                h5_stats_path = f'{QA_PATH}/statistics/data/{pol}/{pwr_phase}'
                stats_group = cslc_h5py_root.require_group(h5_stats_path)

                # build list of QA stat items for pwr_phase
                qa_items = []
                vals = _compute_slc_array_stats(pol_arr, pwr_phase)
                for val_name, val in zip(self.stat_names, vals):
                    desc = f'{val_name} of {pwr_phase} of {pol} geocoded SLC'
                    qa_items.append(Meta(val_name, val, desc))

                # save stats to dict and write to HDF5
                _qa_items_to_h5_and_dict(stats_group, pol_dict[pwr_phase],
                                         qa_items)


    def compute_static_layer_stats(self, cslc_h5py_root, rdr2geo_params):
        '''
        Compute correction stats. Stats written to HDF5 and saved to class dict
        for later JSON output

        Parameters
        ----------
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5
        apply_tropo_corrections: bool
            Whether or not to compute troposphere correction stats
        tropo_delay_type: str
            Type of troposphere delay. Any between 'dry', or 'wet', or
            'wet_dry' for the sum of wet and dry troposphere delays. Only used
            apply_tropo_corrections is true.
        '''
        # path to source group
        static_layer_path = f'{DATA_PATH}'

        # Get the static layer to compute stats
        # Following dict tracks which static layers to generate
        # key: file name of static layer
        # value: bool flag from runconfig that determines if layer is computed
        static_layers_dict = {
            file_name_x: rdr2geo_params.compute_longitude,
            file_name_y: rdr2geo_params.compute_latitude,
            file_name_z: rdr2geo_params.compute_height,
            file_name_local_incidence: rdr2geo_params.compute_local_incidence_angle,
            file_name_los_east: rdr2geo_params.compute_ground_to_sat_east,
            file_name_los_north: rdr2geo_params.compute_ground_to_sat_north
        }
        static_layers = [key for key, val in static_layers_dict.items()
                         if val]

        self.compute_stats_from_float_hdf5_dataset(cslc_h5py_root,
                                                   static_layer_path,
                                                   'static_layers',
                                                   static_layers)


    def compute_correction_stats(self, cslc_h5py_root, apply_tropo_corrections,
                                 tropo_delay_type):
        '''
        Compute correction stats. Stats written to HDF5 and saved to class dict
        for later JSON output

        Parameters
        ----------
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5
        apply_tropo_corrections: bool
            Whether or not to compute troposphere correction stats
        tropo_delay_type: str
            Type of troposphere delay. Any between 'dry', or 'wet', or
            'wet_dry' for the sum of wet and dry troposphere delays. Only used
            apply_tropo_corrections is true.
        '''
        # path to source group
        corrections_src_path = f'{METADATA_PATH}/processing_information/timing_corrections'

        # names of datasets to compute stats for
        corrections = ['bistatic_delay', 'geometry_steering_doppler',
                       'azimuth_fm_rate_mismatch', 'los_ionospheric_delay',
                       'los_solid_earth_tides', 'azimuth_solid_earth_tides']

        # check if tropo corrections need to be computed and saved
        if apply_tropo_corrections:
            for delay_type in ['wet', 'dry']:
                if delay_type in tropo_delay_type:
                    corrections.append(f'{delay_type}_los_troposphere_delay')

        self.compute_stats_from_float_hdf5_dataset(cslc_h5py_root,
                                                   corrections_src_path,
                                                   'timing_corrections', corrections)


    def compute_stats_from_float_hdf5_dataset(self, cslc_h5py_root,
                                              src_group_path, qa_group_name,
                                              qa_item_names):
        '''
        Compute correction stats for float-type, HDF5datasets. Stats written to
        HDF5 and saved to class dict for later JSON output

        Parameters
        ----------
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5
        src_group_path: str
        qa_group_name: str
        qa_item_names: list[str]
        '''
        # init dict to save all QA item stats to
        self.stats_dict[qa_group_name] = {}
        qa_dict = self.stats_dict[qa_group_name]

        # compute stats and write to hdf5
        for qa_item_name in qa_item_names:
            # init dict for current QA item
            qa_dict[qa_item_name] = {}
            qa_item_dict = qa_dict[qa_item_name]

            # get dataset and compute stats according to dtype
            qa_item_path = f'{src_group_path}/{qa_item_name}'
            qa_item_ds = cslc_h5py_root[qa_item_path]

            # compute stats
            stat_obj = isce3.math.StatsFloat32(qa_item_ds[()].astype(np.float32))

            # create HDF5 group for stats of current QA item
            h5_stats_path = f'{QA_PATH}/statistics/{qa_group_name}/{qa_item_name}'
            qa_item_stats_group = cslc_h5py_root.require_group(h5_stats_path)

            # build list of QA stat items
            qa_items = []
            vals = [stat_obj.mean, stat_obj.min, stat_obj.max,
                    stat_obj.sample_stddev]
            for val_name, val in zip(self.stat_names, vals):
                desc = f'{val_name} of {qa_item_name}'
                qa_items.append(Meta(val_name, val, desc))

            # save stats to dict and write to HDF5
            _qa_items_to_h5_and_dict(qa_item_stats_group, qa_item_dict,
                                     qa_items)


    def shadow_pixel_classification(self, cslc_h5py_root):
        '''
        Populate classification of shadow layover pixels

        Parameters
        ----------
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5
        '''

        percent_shadow, percent_layover, percent_combined =\
            self.compute_layover_shadow_pixel_percent(cslc_h5py_root)
        pxl_qa_items = [
            Meta('percent_layover_pixels', percent_layover,
                 'Percentage of output pixels labeled layover'),
            Meta('percent_shadow_pixels', percent_shadow,
                 'Percentage of output pixels labeled shadow'),
            Meta('percent_combined_pixels', percent_combined,
                 'Percentage of output pixels labeled layover and shadow')
        ]

        # create HDF5 group for pixel classification info
        h5_pxl_path = f'{QA_PATH}/pixel_classification'
        pxl_group = cslc_h5py_root.require_group(h5_pxl_path)

        # write items to HDF5 and dict
        _qa_items_to_h5_and_dict(pxl_group, self.pixel_percentage_dict,
                                 pxl_qa_items)


    def percent_land_and_valid_pixels(self, cslc_h5py_root, pol):
        '''
        Populate classification of geocoded pixel types

        Parameters
        ----------
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5

        pol: str
            Polarization of the CSLC layer
        '''
        percent_land_pixels, percent_valid_pixels = \
            self.compute_valid_land_and_pixel_percents(cslc_h5py_root,
                                                       pol)
        pxl_qa_items = [
            Meta('percent_land_pixels', percent_land_pixels,
                 'Percentage of output pixels labeled as land'),
            Meta('percent_valid_pixels', percent_valid_pixels,
                 'Percentage of output pixels are valid')
        ]

        # create HDF5 group for pixel classification info
        h5_pxl_path = f'{QA_PATH}/pixel_classification'
        pxl_group = cslc_h5py_root.require_group(h5_pxl_path)

        # write items to HDF5 and dict
        _qa_items_to_h5_and_dict(pxl_group, self.pixel_percentage_dict,
                                 pxl_qa_items)


    def populate_rfi_dict(self, cslc_h5py_root, bursts):
        '''
        Place holder for populating SAFE RFI information

        Parameters
        ----------
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5
        bursts: list[Sentinel1BurstSlc]
            List of burst SLC object with RFI info
        '''

        for burst in bursts:
            is_rfi_info_available = burst.burst_rfi_info is not None
            rfi_qa_items_pol = [Meta('is_rfi_info_available',
                                is_rfi_info_available,
                                'Whether or not RFI information is available')]

            if is_rfi_info_available:
                # Follow key/values only assigned if RFI info is available
                rfi_info_list = [
                    Meta('rfi_mitigation_performed',
                         burst.burst_rfi_info.rfi_mitigation_performed,
                         ('Activation strategy of RFI mitigation'
                          '["never", "BasedOnNoiseMeas", "always"]')),
                    Meta('rfi_mitigation_domain',
                         burst.burst_rfi_info.rfi_mitigation_domain,
                         'Domain the RFI mitigation step was performed')
                ]
                rfi_qa_items_pol += rfi_info_list

            # create HDF5 group for RFI info for current polarization
            h5_rfi_path = f'{QA_PATH}/rfi_information/{burst.polarization}'
            rfi_group = cslc_h5py_root.require_group(h5_rfi_path)

            # write items to HDF5 and dict
            _qa_items_to_h5_and_dict(rfi_group, self.rfi_dict, rfi_qa_items_pol)

            # Take care of the burst RFI report information
            if not is_rfi_info_available:
                return

            # Alias for readability
            rfi_burst_report = burst.burst_rfi_info.rfi_burst_report

            # Add the metadata of the burst RFI report
            rfi_burst_report_list = [
                Meta('swath',
                     rfi_burst_report['swath'],
                     'Swath of the burst'),
                Meta('azimuth_time',
                     datetime.datetime.strftime(rfi_burst_report['azimuthTime'],
                                                TIME_STR_FMT),
                    'Azimuth time of the burst report'),
                Meta('in_band_out_band_power_ratio',
                          rfi_burst_report['inBandOutBandPowerRatio'],
                          'Ratio between the in-band and out-of-band power of the burst')
            ]

            self.rfi_dict['rfi_burst_report'] = {}
            rfi_burst_report_group = rfi_group.require_group('rfi_burst_report')
            _qa_items_to_h5_and_dict(rfi_burst_report_group,
                                     self.rfi_dict['rfi_burst_report'],
                                     rfi_burst_report_list)

            # Take care of the time domain portion of the burst report
            if 'timeDomainRfiReport' in rfi_burst_report.keys():
                time_domain_report = rfi_burst_report['timeDomainRfiReport']
                burst_time_domain_report_item = [
                    Meta('percentage_affected_lines',
                         time_domain_report['percentageAffectedLines'],
                         'Percentage of level-0 lines affected by RFI.'),
                    Meta('avg_percentage_affected_samples',
                         time_domain_report['avgPercentageAffectedSamples'],
                         'Average percentage of affected level-0 samples in the lines containing RFI'),
                    Meta('max_percentage_affected_samples',
                         time_domain_report['maxPercentageAffectedSamples'],
                         'Maximum percentage of level-0 samples affected by RFI in the same line'),
                ]

                self.rfi_dict['rfi_burst_report']['time_domain_rfi_report'] = {}
                rfi_burst_report_time_domain_group =\
                    rfi_burst_report_group.require_group('time_domain_rfi_report')
                _qa_items_to_h5_and_dict(rfi_burst_report_time_domain_group,
                                         self.rfi_dict['rfi_burst_report']['time_domain_rfi_report'],
                                         burst_time_domain_report_item)

            # Take care of the frequency time domain portion of the burst report
            if 'frequencyDomainRfiBurstReport' in rfi_burst_report.keys():
                freq_domain_report = rfi_burst_report['frequencyDomainRfiBurstReport']
                burst_freq_domain_report_item = [
                    Meta('num_sub_blocks',
                         freq_domain_report['numSubBlocks'],
                         'Number of sub-blocks in the current burst'),
                    Meta('sub_block_size',
                         freq_domain_report['subBlockSize'],
                         'Number of lines in each sub-block'),
                    Meta('percentage_blocks_persistent_rfi',
                         freq_domain_report['percentageBlocksPersistentRfi'],
                         ('Percentage of processing blocks affected by persistent RFI. '
                          'In this case the RFI detection is performed on the mean PSD of '
                          'each processing block')),
                    Meta('max_percentage_bw_affected_persistent_rfi',
                         freq_domain_report['maxPercentageBWAffectedPersistentRfi'],
                         ('Max percentage bandwidth affected by '
                          'persistent RFI in a single processing block.'))
                ]

                self.rfi_dict['rfi_burst_report']['frequency_domain_rfi_report'] = {}
                rfi_burst_report_freq_domain_group = rfi_burst_report_group.require_group('frequency_domain_rfi_report')
                _qa_items_to_h5_and_dict(rfi_burst_report_freq_domain_group,
                                        self.rfi_dict['rfi_burst_report']['frequency_domain_rfi_report'],
                                        burst_freq_domain_report_item)

                # Take care of isolated RFI report inside frequency burst RFI report
                isolated_rfi_report = freq_domain_report['isolatedRfiReport']
                isolated_report_item = [
                    Meta('percentage_affected_lines',
                         isolated_rfi_report['percentageAffectedLines'],
                         'Percentage of level-0 lines affected by isolated RFI'),
                    Meta('max_percentage_affected_bw',
                         isolated_rfi_report['maxPercentageAffectedBW'],
                         'Max. percentage of bandwidth affected by isolated RFI in a single line')
                ]

                self.rfi_dict['rfi_burst_report']['time_domain_rfi_report']['isolated_rfi_report'] = {}
                isolated_rfi_report_group = rfi_burst_report_freq_domain_group.require_group('isolated_rfi_report')
                _qa_items_to_h5_and_dict(isolated_rfi_report_group,
                                        self.rfi_dict['rfi_burst_report']['time_domain_rfi_report']['isolated_rfi_report'],
                                        isolated_report_item)


    def set_orbit_type(self, cfg, cslc_h5py_root):
        '''
        Populate QA orbit information

        Parameters
        ----------
        cfg: dict
            Runconfig dict containing orbit path
        cslc_h5py_root: h5py.File
            Root of CSLC HDF5
        '''
        orbit_file_path = Path(cfg.orbit_path[0]).name
        if 'RESORB' in orbit_file_path:
            orbit_type = 'restituted orbit file'
        if 'POEORB' in orbit_file_path:
            orbit_type = 'precise orbit file'
        orbit_qa_items = [
            Meta('orbit_type', orbit_type,
                 'Type of orbit file used for processing. '
                 'RESORB: restituted orbit ephemeris or POEORB: precise orbit ephemeris')
        ]

        # create HDF5 group for orbit info
        h5_orbit_path = f'{QA_PATH}/orbit_information'
        orbit_group = cslc_h5py_root.require_group(h5_orbit_path)

        # write to HDF5 group orbit info
        _qa_items_to_h5_and_dict(orbit_group, self.orbit_dict, orbit_qa_items)


    def write_qa_dicts_to_json(self, file_path):
        '''
        Write computed stats in dict to JSON file

        Parameters
        ----------
        file_path: str
            JSON file to write stats to
        '''
        # combine all the dicts into one for output
        output_dict = {
            'raster_statistics': self.stats_dict,
            'pixel_classification_percentatges': self.pixel_percentage_dict,
            'rfi_information': self.rfi_dict,
            'orbit_information': self.orbit_dict}

        # write combined dict to JSON
        with open(file_path, 'w') as f:
            json.dump(output_dict, f, indent=4)


    def compute_valid_land_and_pixel_percents(self, cslc_h5py_root, pol):
        '''
        Compute the percentage of valid pixels on land area

        Parameters
        ----------
        cslc_h5py_path: h5py.File
            Root of the CSLC-S1 HDF5 product
        pol: str
            Polarization of CSLC layer to
            compute the valid pixel area

        Returns
        -------
        percent_valid_land_px: float
            Percentage of valid pixels on land
            in the geocoded burst area
        percent_valid_px: float
            Percentage of invalid pixels
            in the geocoded burst area
        '''
        # extract the geogrid information
        epsg_cslc = int(cslc_h5py_root[f'{DATA_PATH}/projection'][()])

        x_spacing = float(cslc_h5py_root[f'{DATA_PATH}/x_spacing'][()])
        y_spacing = float(cslc_h5py_root[f'{DATA_PATH}/y_spacing'][()])

        x0 = list(cslc_h5py_root[f'{DATA_PATH}/x_coordinates'][()])[0] - x_spacing / 2
        y0 = list(cslc_h5py_root[f'{DATA_PATH}/y_coordinates'][()])[0] - y_spacing / 2

        cslc_array = np.array(cslc_h5py_root[f'{DATA_PATH}/{pol}'])

        height_cslc, width_cslc = cslc_array.shape

        mask_land = _get_land_mask(epsg_cslc,
                                   (x0, x_spacing, 0, y0, 0, y_spacing),
                                   (height_cslc, width_cslc))

        mask_geocoded_burst = _get_valid_pixel_mask(cslc_array)

        mask_valid_inside_burst = mask_geocoded_burst & ~np.isnan(cslc_array)

        mask_valid_land_pixel = mask_geocoded_burst & mask_land

        n_unmasked_pxls = mask_geocoded_burst.sum()
        percent_valid_land_px = mask_valid_land_pixel.sum() / n_unmasked_pxls * 100
        percent_valid_px = mask_valid_inside_burst.sum() / n_unmasked_pxls * 100

        return percent_valid_land_px, percent_valid_px


    def compute_layover_shadow_pixel_percent(self, cslc_h5py_root):
        '''
        Compute the percentage of layover, shadow, and
        layover/shadow pixels in the geocoded burst area

        Parameters
        ----------
        cslc_h5py_path: h5py.File
            Root of the CSLC-S1 HDF5 product

        Returns
        -------
        percent_shadow: float
            Percentage of the shadow pixels
            in the geocoded burst area
        percent_layover: float
            Percentage of the layover pixels
            in the geocoded burst area
        percent_combined: float
            Percentage of the shadow and layover pixels
            in the geocoded burst area
        '''
        layover_shadow_mask_array = cslc_h5py_root[f'{DATA_PATH}/layover_shadow_mask'][()]

        mask_geocoded_burst = layover_shadow_mask_array != 127

        n_unmasked_pxls = mask_geocoded_burst.sum()

        mask_shadow_inside_burst = mask_geocoded_burst & (layover_shadow_mask_array == 1)
        percent_shadow = mask_shadow_inside_burst.sum() / n_unmasked_pxls * 100

        mask_layover_inside_burst = mask_geocoded_burst & (layover_shadow_mask_array == 2)
        percent_layover = mask_layover_inside_burst.sum() / n_unmasked_pxls * 100

        mask_combined_inside_burst = mask_geocoded_burst & (layover_shadow_mask_array == 3)
        percent_combined = mask_combined_inside_burst.sum() / n_unmasked_pxls * 100

        return percent_shadow, percent_layover, percent_combined

compute_CSLC_raster_stats(cslc_h5py_root, bursts)

Compute CSLC raster stats. Stats written to HDF5 and saved to class dict for later JSON output

Parameters:

Name Type Description Default
cslc_h5py_root

Root of CSLC HDF5

 required
bursts

Bursts whose geocoded raster stats are to be computed

 required
Source code in src/compass/s1_cslc_qa.py 
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def compute_CSLC_raster_stats(self, cslc_h5py_root, bursts):
    '''
    Compute CSLC raster stats. Stats written to HDF5 and saved to class
    dict for later JSON output

    Parameters
    ----------
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5
    bursts: list
        Bursts whose geocoded raster stats are to be computed
    '''
    for b in bursts:
        pol = b.polarization

        # get dataset and compute stats according to dtype
        pol_path = f'{DATA_PATH}/{pol}'
        pol_arr = cslc_h5py_root[pol_path][()]

        # create dict for current polarization
        self.stats_dict[pol] = {}
        pol_dict = self.stats_dict[pol]

        # compute power or phase then write stats to HDF5 for CSLC
        for pwr_phase in ['power', 'phase']:
            # create dict to store real/imaginary stat items
            pol_dict[pwr_phase] = {}

            # create HDF5 group for power or phase stats of current
            # polarization
            h5_stats_path = f'{QA_PATH}/statistics/data/{pol}/{pwr_phase}'
            stats_group = cslc_h5py_root.require_group(h5_stats_path)

            # build list of QA stat items for pwr_phase
            qa_items = []
            vals = _compute_slc_array_stats(pol_arr, pwr_phase)
            for val_name, val in zip(self.stat_names, vals):
                desc = f'{val_name} of {pwr_phase} of {pol} geocoded SLC'
                qa_items.append(Meta(val_name, val, desc))

            # save stats to dict and write to HDF5
            _qa_items_to_h5_and_dict(stats_group, pol_dict[pwr_phase],
                                     qa_items)

compute_correction_stats(cslc_h5py_root, apply_tropo_corrections, tropo_delay_type)

Compute correction stats. Stats written to HDF5 and saved to class dict for later JSON output

Parameters:

Name Type Description Default
cslc_h5py_root

Root of CSLC HDF5

 required
apply_tropo_corrections

Whether or not to compute troposphere correction stats

 required
tropo_delay_type

Type of troposphere delay. Any between 'dry', or 'wet', or 'wet_dry' for the sum of wet and dry troposphere delays. Only used apply_tropo_corrections is true.

 required
Source code in src/compass/s1_cslc_qa.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def compute_correction_stats(self, cslc_h5py_root, apply_tropo_corrections,
                             tropo_delay_type):
    '''
    Compute correction stats. Stats written to HDF5 and saved to class dict
    for later JSON output

    Parameters
    ----------
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5
    apply_tropo_corrections: bool
        Whether or not to compute troposphere correction stats
    tropo_delay_type: str
        Type of troposphere delay. Any between 'dry', or 'wet', or
        'wet_dry' for the sum of wet and dry troposphere delays. Only used
        apply_tropo_corrections is true.
    '''
    # path to source group
    corrections_src_path = f'{METADATA_PATH}/processing_information/timing_corrections'

    # names of datasets to compute stats for
    corrections = ['bistatic_delay', 'geometry_steering_doppler',
                   'azimuth_fm_rate_mismatch', 'los_ionospheric_delay',
                   'los_solid_earth_tides', 'azimuth_solid_earth_tides']

    # check if tropo corrections need to be computed and saved
    if apply_tropo_corrections:
        for delay_type in ['wet', 'dry']:
            if delay_type in tropo_delay_type:
                corrections.append(f'{delay_type}_los_troposphere_delay')

    self.compute_stats_from_float_hdf5_dataset(cslc_h5py_root,
                                               corrections_src_path,
                                               'timing_corrections', corrections)

compute_layover_shadow_pixel_percent(cslc_h5py_root)

Compute the percentage of layover, shadow, and layover/shadow pixels in the geocoded burst area

Parameters:

Name Type Description Default
cslc_h5py_path

Root of the CSLC-S1 HDF5 product

 required

Returns:

Name Type Description
percent_shadow float

Percentage of the shadow pixels in the geocoded burst area
percent_layover float

Percentage of the layover pixels in the geocoded burst area
percent_combined float

Percentage of the shadow and layover pixels in the geocoded burst area
Source code in src/compass/s1_cslc_qa.py 
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
def compute_layover_shadow_pixel_percent(self, cslc_h5py_root):
    '''
    Compute the percentage of layover, shadow, and
    layover/shadow pixels in the geocoded burst area

    Parameters
    ----------
    cslc_h5py_path: h5py.File
        Root of the CSLC-S1 HDF5 product

    Returns
    -------
    percent_shadow: float
        Percentage of the shadow pixels
        in the geocoded burst area
    percent_layover: float
        Percentage of the layover pixels
        in the geocoded burst area
    percent_combined: float
        Percentage of the shadow and layover pixels
        in the geocoded burst area
    '''
    layover_shadow_mask_array = cslc_h5py_root[f'{DATA_PATH}/layover_shadow_mask'][()]

    mask_geocoded_burst = layover_shadow_mask_array != 127

    n_unmasked_pxls = mask_geocoded_burst.sum()

    mask_shadow_inside_burst = mask_geocoded_burst & (layover_shadow_mask_array == 1)
    percent_shadow = mask_shadow_inside_burst.sum() / n_unmasked_pxls * 100

    mask_layover_inside_burst = mask_geocoded_burst & (layover_shadow_mask_array == 2)
    percent_layover = mask_layover_inside_burst.sum() / n_unmasked_pxls * 100

    mask_combined_inside_burst = mask_geocoded_burst & (layover_shadow_mask_array == 3)
    percent_combined = mask_combined_inside_burst.sum() / n_unmasked_pxls * 100

    return percent_shadow, percent_layover, percent_combined

compute_static_layer_stats(cslc_h5py_root, rdr2geo_params)

Compute correction stats. Stats written to HDF5 and saved to class dict for later JSON output

Parameters:

Name Type Description Default
cslc_h5py_root

Root of CSLC HDF5

 required
apply_tropo_corrections

Whether or not to compute troposphere correction stats

 required
tropo_delay_type

Type of troposphere delay. Any between 'dry', or 'wet', or 'wet_dry' for the sum of wet and dry troposphere delays. Only used apply_tropo_corrections is true.

 required
Source code in src/compass/s1_cslc_qa.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
def compute_static_layer_stats(self, cslc_h5py_root, rdr2geo_params):
    '''
    Compute correction stats. Stats written to HDF5 and saved to class dict
    for later JSON output

    Parameters
    ----------
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5
    apply_tropo_corrections: bool
        Whether or not to compute troposphere correction stats
    tropo_delay_type: str
        Type of troposphere delay. Any between 'dry', or 'wet', or
        'wet_dry' for the sum of wet and dry troposphere delays. Only used
        apply_tropo_corrections is true.
    '''
    # path to source group
    static_layer_path = f'{DATA_PATH}'

    # Get the static layer to compute stats
    # Following dict tracks which static layers to generate
    # key: file name of static layer
    # value: bool flag from runconfig that determines if layer is computed
    static_layers_dict = {
        file_name_x: rdr2geo_params.compute_longitude,
        file_name_y: rdr2geo_params.compute_latitude,
        file_name_z: rdr2geo_params.compute_height,
        file_name_local_incidence: rdr2geo_params.compute_local_incidence_angle,
        file_name_los_east: rdr2geo_params.compute_ground_to_sat_east,
        file_name_los_north: rdr2geo_params.compute_ground_to_sat_north
    }
    static_layers = [key for key, val in static_layers_dict.items()
                     if val]

    self.compute_stats_from_float_hdf5_dataset(cslc_h5py_root,
                                               static_layer_path,
                                               'static_layers',
                                               static_layers)

compute_stats_from_float_hdf5_dataset(cslc_h5py_root, src_group_path, qa_group_name, qa_item_names)

Compute correction stats for float-type, HDF5datasets. Stats written to HDF5 and saved to class dict for later JSON output

Parameters:

Name Type Description Default
cslc_h5py_root

Root of CSLC HDF5

 required
src_group_path
required
qa_group_name
required
qa_item_names
required
Source code in src/compass/s1_cslc_qa.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def compute_stats_from_float_hdf5_dataset(self, cslc_h5py_root,
                                          src_group_path, qa_group_name,
                                          qa_item_names):
    '''
    Compute correction stats for float-type, HDF5datasets. Stats written to
    HDF5 and saved to class dict for later JSON output

    Parameters
    ----------
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5
    src_group_path: str
    qa_group_name: str
    qa_item_names: list[str]
    '''
    # init dict to save all QA item stats to
    self.stats_dict[qa_group_name] = {}
    qa_dict = self.stats_dict[qa_group_name]

    # compute stats and write to hdf5
    for qa_item_name in qa_item_names:
        # init dict for current QA item
        qa_dict[qa_item_name] = {}
        qa_item_dict = qa_dict[qa_item_name]

        # get dataset and compute stats according to dtype
        qa_item_path = f'{src_group_path}/{qa_item_name}'
        qa_item_ds = cslc_h5py_root[qa_item_path]

        # compute stats
        stat_obj = isce3.math.StatsFloat32(qa_item_ds[()].astype(np.float32))

        # create HDF5 group for stats of current QA item
        h5_stats_path = f'{QA_PATH}/statistics/{qa_group_name}/{qa_item_name}'
        qa_item_stats_group = cslc_h5py_root.require_group(h5_stats_path)

        # build list of QA stat items
        qa_items = []
        vals = [stat_obj.mean, stat_obj.min, stat_obj.max,
                stat_obj.sample_stddev]
        for val_name, val in zip(self.stat_names, vals):
            desc = f'{val_name} of {qa_item_name}'
            qa_items.append(Meta(val_name, val, desc))

        # save stats to dict and write to HDF5
        _qa_items_to_h5_and_dict(qa_item_stats_group, qa_item_dict,
                                 qa_items)

compute_valid_land_and_pixel_percents(cslc_h5py_root, pol)

Compute the percentage of valid pixels on land area

Parameters:

Name Type Description Default
cslc_h5py_path

Root of the CSLC-S1 HDF5 product

 required
pol

Polarization of CSLC layer to compute the valid pixel area

 required

Returns:

Name Type Description
percent_valid_land_px float

Percentage of valid pixels on land in the geocoded burst area
percent_valid_px float

Percentage of invalid pixels in the geocoded burst area
Source code in src/compass/s1_cslc_qa.py 
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
def compute_valid_land_and_pixel_percents(self, cslc_h5py_root, pol):
    '''
    Compute the percentage of valid pixels on land area

    Parameters
    ----------
    cslc_h5py_path: h5py.File
        Root of the CSLC-S1 HDF5 product
    pol: str
        Polarization of CSLC layer to
        compute the valid pixel area

    Returns
    -------
    percent_valid_land_px: float
        Percentage of valid pixels on land
        in the geocoded burst area
    percent_valid_px: float
        Percentage of invalid pixels
        in the geocoded burst area
    '''
    # extract the geogrid information
    epsg_cslc = int(cslc_h5py_root[f'{DATA_PATH}/projection'][()])

    x_spacing = float(cslc_h5py_root[f'{DATA_PATH}/x_spacing'][()])
    y_spacing = float(cslc_h5py_root[f'{DATA_PATH}/y_spacing'][()])

    x0 = list(cslc_h5py_root[f'{DATA_PATH}/x_coordinates'][()])[0] - x_spacing / 2
    y0 = list(cslc_h5py_root[f'{DATA_PATH}/y_coordinates'][()])[0] - y_spacing / 2

    cslc_array = np.array(cslc_h5py_root[f'{DATA_PATH}/{pol}'])

    height_cslc, width_cslc = cslc_array.shape

    mask_land = _get_land_mask(epsg_cslc,
                               (x0, x_spacing, 0, y0, 0, y_spacing),
                               (height_cslc, width_cslc))

    mask_geocoded_burst = _get_valid_pixel_mask(cslc_array)

    mask_valid_inside_burst = mask_geocoded_burst & ~np.isnan(cslc_array)

    mask_valid_land_pixel = mask_geocoded_burst & mask_land

    n_unmasked_pxls = mask_geocoded_burst.sum()
    percent_valid_land_px = mask_valid_land_pixel.sum() / n_unmasked_pxls * 100
    percent_valid_px = mask_valid_inside_burst.sum() / n_unmasked_pxls * 100

    return percent_valid_land_px, percent_valid_px

percent_land_and_valid_pixels(cslc_h5py_root, pol)

Populate classification of geocoded pixel types

Parameters:

Name Type Description Default
cslc_h5py_root

Root of CSLC HDF5

 required
pol

Polarization of the CSLC layer

 required
Source code in src/compass/s1_cslc_qa.py 
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
def percent_land_and_valid_pixels(self, cslc_h5py_root, pol):
    '''
    Populate classification of geocoded pixel types

    Parameters
    ----------
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5

    pol: str
        Polarization of the CSLC layer
    '''
    percent_land_pixels, percent_valid_pixels = \
        self.compute_valid_land_and_pixel_percents(cslc_h5py_root,
                                                   pol)
    pxl_qa_items = [
        Meta('percent_land_pixels', percent_land_pixels,
             'Percentage of output pixels labeled as land'),
        Meta('percent_valid_pixels', percent_valid_pixels,
             'Percentage of output pixels are valid')
    ]

    # create HDF5 group for pixel classification info
    h5_pxl_path = f'{QA_PATH}/pixel_classification'
    pxl_group = cslc_h5py_root.require_group(h5_pxl_path)

    # write items to HDF5 and dict
    _qa_items_to_h5_and_dict(pxl_group, self.pixel_percentage_dict,
                             pxl_qa_items)

populate_rfi_dict(cslc_h5py_root, bursts)

Place holder for populating SAFE RFI information

Parameters:

Name Type Description Default
cslc_h5py_root

Root of CSLC HDF5

 required
bursts

List of burst SLC object with RFI info

 required
Source code in src/compass/s1_cslc_qa.py 
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
def populate_rfi_dict(self, cslc_h5py_root, bursts):
    '''
    Place holder for populating SAFE RFI information

    Parameters
    ----------
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5
    bursts: list[Sentinel1BurstSlc]
        List of burst SLC object with RFI info
    '''

    for burst in bursts:
        is_rfi_info_available = burst.burst_rfi_info is not None
        rfi_qa_items_pol = [Meta('is_rfi_info_available',
                            is_rfi_info_available,
                            'Whether or not RFI information is available')]

        if is_rfi_info_available:
            # Follow key/values only assigned if RFI info is available
            rfi_info_list = [
                Meta('rfi_mitigation_performed',
                     burst.burst_rfi_info.rfi_mitigation_performed,
                     ('Activation strategy of RFI mitigation'
                      '["never", "BasedOnNoiseMeas", "always"]')),
                Meta('rfi_mitigation_domain',
                     burst.burst_rfi_info.rfi_mitigation_domain,
                     'Domain the RFI mitigation step was performed')
            ]
            rfi_qa_items_pol += rfi_info_list

        # create HDF5 group for RFI info for current polarization
        h5_rfi_path = f'{QA_PATH}/rfi_information/{burst.polarization}'
        rfi_group = cslc_h5py_root.require_group(h5_rfi_path)

        # write items to HDF5 and dict
        _qa_items_to_h5_and_dict(rfi_group, self.rfi_dict, rfi_qa_items_pol)

        # Take care of the burst RFI report information
        if not is_rfi_info_available:
            return

        # Alias for readability
        rfi_burst_report = burst.burst_rfi_info.rfi_burst_report

        # Add the metadata of the burst RFI report
        rfi_burst_report_list = [
            Meta('swath',
                 rfi_burst_report['swath'],
                 'Swath of the burst'),
            Meta('azimuth_time',
                 datetime.datetime.strftime(rfi_burst_report['azimuthTime'],
                                            TIME_STR_FMT),
                'Azimuth time of the burst report'),
            Meta('in_band_out_band_power_ratio',
                      rfi_burst_report['inBandOutBandPowerRatio'],
                      'Ratio between the in-band and out-of-band power of the burst')
        ]

        self.rfi_dict['rfi_burst_report'] = {}
        rfi_burst_report_group = rfi_group.require_group('rfi_burst_report')
        _qa_items_to_h5_and_dict(rfi_burst_report_group,
                                 self.rfi_dict['rfi_burst_report'],
                                 rfi_burst_report_list)

        # Take care of the time domain portion of the burst report
        if 'timeDomainRfiReport' in rfi_burst_report.keys():
            time_domain_report = rfi_burst_report['timeDomainRfiReport']
            burst_time_domain_report_item = [
                Meta('percentage_affected_lines',
                     time_domain_report['percentageAffectedLines'],
                     'Percentage of level-0 lines affected by RFI.'),
                Meta('avg_percentage_affected_samples',
                     time_domain_report['avgPercentageAffectedSamples'],
                     'Average percentage of affected level-0 samples in the lines containing RFI'),
                Meta('max_percentage_affected_samples',
                     time_domain_report['maxPercentageAffectedSamples'],
                     'Maximum percentage of level-0 samples affected by RFI in the same line'),
            ]

            self.rfi_dict['rfi_burst_report']['time_domain_rfi_report'] = {}
            rfi_burst_report_time_domain_group =\
                rfi_burst_report_group.require_group('time_domain_rfi_report')
            _qa_items_to_h5_and_dict(rfi_burst_report_time_domain_group,
                                     self.rfi_dict['rfi_burst_report']['time_domain_rfi_report'],
                                     burst_time_domain_report_item)

        # Take care of the frequency time domain portion of the burst report
        if 'frequencyDomainRfiBurstReport' in rfi_burst_report.keys():
            freq_domain_report = rfi_burst_report['frequencyDomainRfiBurstReport']
            burst_freq_domain_report_item = [
                Meta('num_sub_blocks',
                     freq_domain_report['numSubBlocks'],
                     'Number of sub-blocks in the current burst'),
                Meta('sub_block_size',
                     freq_domain_report['subBlockSize'],
                     'Number of lines in each sub-block'),
                Meta('percentage_blocks_persistent_rfi',
                     freq_domain_report['percentageBlocksPersistentRfi'],
                     ('Percentage of processing blocks affected by persistent RFI. '
                      'In this case the RFI detection is performed on the mean PSD of '
                      'each processing block')),
                Meta('max_percentage_bw_affected_persistent_rfi',
                     freq_domain_report['maxPercentageBWAffectedPersistentRfi'],
                     ('Max percentage bandwidth affected by '
                      'persistent RFI in a single processing block.'))
            ]

            self.rfi_dict['rfi_burst_report']['frequency_domain_rfi_report'] = {}
            rfi_burst_report_freq_domain_group = rfi_burst_report_group.require_group('frequency_domain_rfi_report')
            _qa_items_to_h5_and_dict(rfi_burst_report_freq_domain_group,
                                    self.rfi_dict['rfi_burst_report']['frequency_domain_rfi_report'],
                                    burst_freq_domain_report_item)

            # Take care of isolated RFI report inside frequency burst RFI report
            isolated_rfi_report = freq_domain_report['isolatedRfiReport']
            isolated_report_item = [
                Meta('percentage_affected_lines',
                     isolated_rfi_report['percentageAffectedLines'],
                     'Percentage of level-0 lines affected by isolated RFI'),
                Meta('max_percentage_affected_bw',
                     isolated_rfi_report['maxPercentageAffectedBW'],
                     'Max. percentage of bandwidth affected by isolated RFI in a single line')
            ]

            self.rfi_dict['rfi_burst_report']['time_domain_rfi_report']['isolated_rfi_report'] = {}
            isolated_rfi_report_group = rfi_burst_report_freq_domain_group.require_group('isolated_rfi_report')
            _qa_items_to_h5_and_dict(isolated_rfi_report_group,
                                    self.rfi_dict['rfi_burst_report']['time_domain_rfi_report']['isolated_rfi_report'],
                                    isolated_report_item)

set_orbit_type(cfg, cslc_h5py_root)

Populate QA orbit information

Parameters:

Name Type Description Default
cfg

Runconfig dict containing orbit path

 required
cslc_h5py_root

Root of CSLC HDF5

 required
Source code in src/compass/s1_cslc_qa.py 
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
def set_orbit_type(self, cfg, cslc_h5py_root):
    '''
    Populate QA orbit information

    Parameters
    ----------
    cfg: dict
        Runconfig dict containing orbit path
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5
    '''
    orbit_file_path = Path(cfg.orbit_path[0]).name
    if 'RESORB' in orbit_file_path:
        orbit_type = 'restituted orbit file'
    if 'POEORB' in orbit_file_path:
        orbit_type = 'precise orbit file'
    orbit_qa_items = [
        Meta('orbit_type', orbit_type,
             'Type of orbit file used for processing. '
             'RESORB: restituted orbit ephemeris or POEORB: precise orbit ephemeris')
    ]

    # create HDF5 group for orbit info
    h5_orbit_path = f'{QA_PATH}/orbit_information'
    orbit_group = cslc_h5py_root.require_group(h5_orbit_path)

    # write to HDF5 group orbit info
    _qa_items_to_h5_and_dict(orbit_group, self.orbit_dict, orbit_qa_items)

shadow_pixel_classification(cslc_h5py_root)

Populate classification of shadow layover pixels

Parameters:

Name Type Description Default
cslc_h5py_root

Root of CSLC HDF5

 required
Source code in src/compass/s1_cslc_qa.py 
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def shadow_pixel_classification(self, cslc_h5py_root):
    '''
    Populate classification of shadow layover pixels

    Parameters
    ----------
    cslc_h5py_root: h5py.File
        Root of CSLC HDF5
    '''

    percent_shadow, percent_layover, percent_combined =\
        self.compute_layover_shadow_pixel_percent(cslc_h5py_root)
    pxl_qa_items = [
        Meta('percent_layover_pixels', percent_layover,
             'Percentage of output pixels labeled layover'),
        Meta('percent_shadow_pixels', percent_shadow,
             'Percentage of output pixels labeled shadow'),
        Meta('percent_combined_pixels', percent_combined,
             'Percentage of output pixels labeled layover and shadow')
    ]

    # create HDF5 group for pixel classification info
    h5_pxl_path = f'{QA_PATH}/pixel_classification'
    pxl_group = cslc_h5py_root.require_group(h5_pxl_path)

    # write items to HDF5 and dict
    _qa_items_to_h5_and_dict(pxl_group, self.pixel_percentage_dict,
                             pxl_qa_items)

write_qa_dicts_to_json(file_path)

Write computed stats in dict to JSON file

Parameters:

Name Type Description Default
file_path

JSON file to write stats to

 required
Source code in src/compass/s1_cslc_qa.py 
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
def write_qa_dicts_to_json(self, file_path):
    '''
    Write computed stats in dict to JSON file

    Parameters
    ----------
    file_path: str
        JSON file to write stats to
    '''
    # combine all the dicts into one for output
    output_dict = {
        'raster_statistics': self.stats_dict,
        'pixel_classification_percentatges': self.pixel_percentage_dict,
        'rfi_information': self.rfi_dict,
        'orbit_information': self.orbit_dict}

    # write combined dict to JSON
    with open(file_path, 'w') as f:
        json.dump(output_dict, f, indent=4)

value_description_dict(val, desc)

Convenience function that returns dict with description and value

Source code in src/compass/s1_cslc_qa.py 
39
40
41
42
43
def value_description_dict(val, desc):
    '''
    Convenience function that returns dict with description and value
    '''
    return {'value': val, 'description': desc}