Merge pull request #15 from vincelhx/main

added a dataset in the datatree for antennaPattern

Author: agrouaze

safe_s1/metadata.py

@@ -100,6 +100,8 @@ def __init__(self, name, backend_kwargs=None):
                 'calibration_luts': self.get_calibration_luts,
                 'noise_azimuth_raw': self.get_noise_azi_raw,
                 'noise_range_raw': self.get_noise_range_raw,
+                'antenna_pattern':self.antenna_pattern,
+                'swath_merging': self.swath_merging
             }
             self.dt = datatree.DataTree.from_dict(self._dict)
             assert self.dt==self.datatree
@@ -496,6 +498,25 @@ def bursts(self):
                                                                        describe=True)
             return bursts
 
+    @property
+    def antenna_pattern(self):
+        ds =  self.xml_parser.get_compound_var(self.files['annotation'].iloc[0], 'antenna_pattern')
+        ds.attrs['history'] = self.xml_parser.get_compound_var(self.files['annotation'].iloc[0], 'antenna_pattern',
+                                                               describe=True)
+        return ds
+
+    @property
+    def swath_merging(self):
+        if 'GRD' in self.product:
+            
+            ds =  self.xml_parser.get_compound_var(self.files['annotation'].iloc[0], 'swath_merging')
+            ds.attrs['history'] = self.xml_parser.get_compound_var(self.files['annotation'].iloc[0], 'swath_merging',
+                                                               describe=True)
+        else :
+            ds = xr.Dataset()
+            
+        return ds
+
     @property
     def multidataset(self):
         """

safe_s1/sentinel1_xml_mappings.py

@@ -148,8 +148,10 @@ def list_poly_from_list_string_coords(str_coords_list):
         'satellite': (scalar, '//safe:platform/safe:number'),
         'start_date': (date_converter, '//safe:acquisitionPeriod/safe:startTime'),
         'stop_date': (date_converter, '//safe:acquisitionPeriod/safe:stopTime'),
-        'aux_cal': (scalar,
-                    '//metadataSection/metadataObject/metadataWrap/xmlData/safe:processing/safe:resource/safe:processing/safe:resource[@role="AUX_CAL"]/@name'),
+        
+        'aux_cal': (scalar, '//metadataSection/metadataObject/metadataWrap/xmlData/safe:processing/safe:resource/safe:processing/safe:resource[@role="AUX_CAL"]/@name'),
+        'aux_pp1': (scalar, '//metadataSection/metadataObject/metadataWrap/xmlData/safe:processing/safe:resource/safe:processing/safe:resource[@role="AUX_PP1"]/@name'),
+
         'aux_cal_sl2': (scalar,'//metadataSection/metadataObject/metadataWrap/xmlData/safe:processing/safe:resource[@role="AUX_CAL"]/@name'),
         'annotation_files': (
             normpath, '/xfdu:XFDU/dataObjectSection/*[@repID="s1Level1ProductSchema"]/byteStream/fileLocation/@href'),
@@ -302,8 +304,32 @@ def list_poly_from_list_string_coords(str_coords_list):
         'fmrate_azimuthFmRatePolynomial': (
             list_of_float_1D_array_from_string,
             '//product/generalAnnotation/azimuthFmRateList/azimuthFmRate/azimuthFmRatePolynomial'),
-
-    },
+        
+        'ap_azimuthTime': (
+            datetime64_array, '/product/antennaPattern/antennaPatternList/antennaPattern/azimuthTime'),        
+        'ap_roll' : (float_array, '/product/antennaPattern/antennaPatternList/antennaPattern/roll'),
+        'ap_swath' : (lambda x: np.array(x), '/product/antennaPattern/antennaPatternList/antennaPattern/swath'),
+        
+        'ap_elevationAngle': (
+            list_of_float_1D_array_from_string, '/product/antennaPattern/antennaPatternList/antennaPattern/elevationAngle'), 
+        'ap_incidenceAngle': (
+            list_of_float_1D_array_from_string, '/product/antennaPattern/antennaPatternList/antennaPattern/incidenceAngle'), 
+        'ap_slantRangeTime': (
+            list_of_float_1D_array_from_string, '/product/antennaPattern/antennaPatternList/antennaPattern/slantRangeTime'),        
+        'ap_terrainHeight': (
+            float_array, '/product/antennaPattern/antennaPatternList/antennaPattern/terrainHeight'),        
+        'ap_elevationPattern' : (
+            list_of_float_1D_array_from_string, '/product/antennaPattern/antennaPatternList/antennaPattern/elevationPattern'),        
+        
+        'sm_nbPerSwat': (int_array, '/product/swathMerging/swathMergeList/swathMerge/swathBoundsList/@count'),
+        'sm_swath' : (lambda x: np.array(x),     '/product/swathMerging/swathMergeList/swathMerge/swath'),
+        'sm_azimuthTime' : (datetime64_array, '/product/swathMerging/swathMergeList/swathMerge/swathBoundsList/swathBounds/azimuthTime'),      
+        'sm_firstAzimuthLine' : (int_array, '/product/swathMerging/swathMergeList/swathMerge/swathBoundsList/swathBounds/firstAzimuthLine'),      
+        'sm_lastAzimuthLine' : (int_array, '/product/swathMerging/swathMergeList/swathMerge/swathBoundsList/swathBounds/lastAzimuthLine'),      
+        'sm_firstRangeSample' : (int_array, '/product/swathMerging/swathMergeList/swathMerge/swathBoundsList/swathBounds/firstRangeSample'),      
+        'sm_lastRangeSample' : (int_array, '/product/swathMerging/swathMergeList/swathMerge/swathBoundsList/swathBounds/lastRangeSample'),      
+
+    },   
     'xsd': {'all': (str, '/xsd:schema/xsd:complexType/xsd:sequence/xsd:element/xsd:annotation/xsd:documentation'),
             'names': (str, '/xsd:schema/xsd:complexType/xsd:sequence/xsd:element/@name'),
             'sensingtime': (str, '/xsd:schema/xsd:complexType/xsd:sequence/xsd:element/sensingTime')
@@ -694,6 +720,7 @@ def doppler_centroid_estimates(nb_dcestimate,
                                                           'source': xpath_mappings['annotation']['dc_rmserrAboveThres'][
                                                               1]},
                                                       coords={'azimuthTime': dc_azimuth_time})
+    
     return ds
 
 
@@ -716,6 +743,147 @@ def geolocation_grid(line, sample, values):
     values = np.reshape(values, shape)
     return xr.DataArray(values, dims=['line', 'sample'], coords={'line': line, 'sample': sample})
 
+def antenna_pattern(ap_swath,ap_roll,ap_azimuthTime,ap_terrainHeight,ap_elevationAngle,ap_elevationPattern,ap_incidenceAngle,ap_slantRangeTime):
+    """
+
+    Parameters
+    ----------
+    ap_swath
+    ap_roll
+    ap_azimuthTime
+    ap_terrainHeight
+    ap_elevationAngle
+    ap_elevationPattern
+    ap_incidenceAngle
+    ap_slantRangeTime
+    
+    Returns
+    -------
+    xarray.DataSet
+    """   
+    # Fonction to convert string 'EW1' ou 'IW3' as int
+    def convert_to_int(swath):
+        return int(swath[-1])
+    vectorized_convert = np.vectorize(convert_to_int)
+    swathNumber = vectorized_convert(ap_swath)
+
+    dim_azimuthTime = max(np.bincount(swathNumber))
+    dim_slantRangeTime = max(array.shape[0] for array in ap_elevationAngle)
+
+    include_roll = len(ap_roll) != 0
+
+    # Create 2Ds arrays 
+    elevAngle2d = np.full((len(ap_elevationAngle), dim_slantRangeTime), np.nan)  
+    gain2d = np.full((len(ap_elevationPattern), dim_slantRangeTime), np.nan)  
+    slantRangeTime2d = np.full((len(ap_slantRangeTime), dim_slantRangeTime), np.nan)
+    incAngle2d = np.full((len(ap_incidenceAngle), dim_slantRangeTime), np.nan)
+
+    
+    for i in range(len(ap_elevationAngle)):
+        elevAngle2d[i, :ap_elevationAngle[i].shape[0]] = ap_elevationAngle[i]
+
+        if ap_elevationAngle[i].shape[0] != ap_elevationPattern[i].shape[0] :
+            gain2d[i, :ap_elevationAngle[i].shape[0]] = np.sqrt(ap_elevationPattern[i][::2]**2+ap_elevationPattern[i][1::2]**2)
+        else:
+            #logging.warn("antenna pattern is not given in complex values. You probably use an old file\n" + e) 
+            gain2d[i, :ap_elevationAngle[i].shape[0]] = ap_elevationPattern[i]
+
+        slantRangeTime2d[i, :ap_slantRangeTime[i].shape[0]] = ap_slantRangeTime[i]
+        incAngle2d[i, :ap_incidenceAngle[i].shape[0]] = ap_incidenceAngle[i]
+
+
+    swath_number_2d = np.full((len(np.unique(swathNumber)), dim_azimuthTime), np.nan)
+    roll_angle_2d = np.full((len(np.unique(swathNumber)), dim_azimuthTime), np.nan)
+    azimuthTime_2d = np.full((len(np.unique(swathNumber)), dim_azimuthTime), np.nan)
+    terrainHeight_2d = np.full((len(np.unique(swathNumber)), dim_azimuthTime), np.nan)
+
+    slantRangeTime_2d = np.full((len(np.unique(swathNumber)), dim_slantRangeTime), np.nan)
+
+    elevationAngle_3d = np.full((len(np.unique(swathNumber)), dim_azimuthTime, dim_slantRangeTime), np.nan)
+    incidenceAngle_3d = np.full((len(np.unique(swathNumber)), dim_azimuthTime, dim_slantRangeTime), np.nan)
+    gain3d = np.full((len(np.unique(swathNumber)), dim_azimuthTime, dim_slantRangeTime), np.nan)
+
+    
+    for i, swath_number in enumerate(np.unique(swathNumber)):
+        length_dim0 = len(ap_azimuthTime[swathNumber == swath_number])
+        swath_number_2d[i, :length_dim0] = swathNumber[swathNumber == swath_number]
+        azimuthTime_2d[i, :length_dim0] = ap_azimuthTime[swathNumber == swath_number]
+        terrainHeight_2d[i, :length_dim0] = ap_terrainHeight[swathNumber == swath_number]
+        slantRangeTime_2d[i, :] = slantRangeTime2d[i, :]
+
+        if include_roll:
+            roll_angle_2d[i, :length_dim0] = ap_roll[swathNumber == swath_number]        
+
+        for j in range(0, dim_slantRangeTime):
+            elevationAngle_3d[i,:length_dim0,j]=elevAngle2d[swathNumber == swath_number,j]
+            incidenceAngle_3d[i,:length_dim0,j]=incAngle2d[swathNumber == swath_number,j]
+            gain3d[i,:length_dim0,j]=gain2d[swathNumber == swath_number,j]
+        
+    azimuthTime_2d = azimuthTime_2d.astype('datetime64[ns]')
+
+    # return a Dataset
+    ds = xr.Dataset({
+        'slantRangeTime' : (['swath_nb', 'dim_slantRangeTime'], slantRangeTime_2d),
+        'swath' : (['swath_nb', 'dim_azimuthTime'], swath_number_2d),
+        'roll' : (['swath_nb', 'dim_azimuthTime'], roll_angle_2d),
+        'azimuthTime' : (['swath_nb', 'dim_azimuthTime'], azimuthTime_2d),
+        'terrainHeight' : (['swath_nb', 'dim_azimuthTime'], terrainHeight_2d),
+        'elevationAngle' : (['swath_nb', 'dim_azimuthTime','dim_slantRangeTime'],elevationAngle_3d),
+        'incidenceAngle' : (['swath_nb', 'dim_azimuthTime','dim_slantRangeTime'],incidenceAngle_3d),
+        'gain' : (['swath_nb', 'dim_azimuthTime','dim_slantRangeTime'],gain3d),
+        },    
+        coords={'swath_nb': np.unique(swathNumber)}
+    )
+    ds.attrs["dim_azimuthTime"] = "max dimension of azimuthTime for a swath"
+    ds.attrs["dim_slantRangeTime"] = "max dimension of slantRangeTime for a swath"
+    ds.attrs["comment"] = "The antenna pattern data set record contains a list of vectors of the \
+                           antenna elevation pattern values that have been updated along track\
+                           and used to correct the radiometry during image processing."
+    ds.attrs["example"] = "for example, if swath Y is smaller than swath X, user has to remove nan to get the dims of the swath"
+    ds.attrs["source"] = "Sentinel-1 Product Specification"
+
+    return ds 
+
+def swath_merging(sm_swath,sm_nbPerSwat,sm_azimuthTime,sm_firstAzimuthLine,sm_lastAzimuthLine,sm_firstRangeSample,sm_lastRangeSample):
+    """
+
+    Parameters
+    ----------
+    sm_swath
+    sm_nbPerSwat
+    sm_azimuthTime
+    sm_firstAzimuthLine
+    sm_lastAzimuthLine
+    sm_firstRangeSample
+    sm_lastRangeSample
+    
+    Returns
+    -------
+    xarray.DataSet
+    """   
+    # Fonction to convert string 'EW1' ou 'IW3' as int
+    def convert_to_int(swath):
+        return int(swath[-1])
+    vectorized_convert = np.vectorize(convert_to_int)
+    repeated_swaths = np.repeat(sm_swath, sm_nbPerSwat)
+    swathNumber = vectorized_convert(repeated_swaths)
+    
+    ds = xr.Dataset({
+        'swaths' : (['dim_azimuthTime'], swathNumber),
+        'azimuthTime' : (['dim_azimuthTime'], sm_azimuthTime),
+        'firstAzimuthLine' : (['dim_azimuthTime'], sm_firstAzimuthLine),
+        'lastAzimuthLine' : (['dim_azimuthTime'], sm_lastAzimuthLine),
+        'firstRangeSample' : (['dim_azimuthTime'], sm_firstRangeSample),
+        'lastRangeSample' : (['dim_azimuthTime'], sm_lastRangeSample),
+        },    
+    )
+    ds.attrs["comment"] = "The swath merging data set record contains information about how \
+                           multiple swaths were stitched together to form one large contiguous \
+                           swath. This data set record only applies to IW and EW GRD \
+                           products"
+    ds.attrs["source"] = "Sentinel-1 Product Specification"
+
+    return ds
 
 # dict of compounds variables.
 # compounds variables are variables composed of several variables.
@@ -732,7 +900,8 @@ def geolocation_grid(line, sample, values):
         'start_date': 'manifest.start_date',
         'stop_date': 'manifest.stop_date',
         'footprints': 'manifest.footprints',
-        'aux_cal': 'manifest.aux_cal'
+        'aux_cal': 'manifest.aux_cal',
+        'aux_pp1': 'manifest.aux_pp1'
     },
     'safe_attributes_sl2': {
         'ipf_version': 'manifest.ipf_version',
@@ -856,4 +1025,17 @@ def geolocation_grid(line, sample, values):
 
                  ),
     },
+    'antenna_pattern': {
+        'func': antenna_pattern,
+        'args': ('annotation.ap_swath','annotation.ap_roll','annotation.ap_azimuthTime','annotation.ap_terrainHeight',
+                 'annotation.ap_elevationAngle','annotation.ap_elevationPattern','annotation.ap_incidenceAngle',
+                 'annotation.ap_slantRangeTime'
+        )
+    },
+    'swath_merging': {
+        'func': swath_merging,
+        'args': ('annotation.sm_swath','annotation.sm_nbPerSwat','annotation.sm_azimuthTime','annotation.sm_firstAzimuthLine',
+                 'annotation.sm_lastAzimuthLine','annotation.sm_firstRangeSample','annotation.sm_lastRangeSample'
+        )
+    },
 }