SAR-ARD · johntruckenbrodt · Aug 25, 2023 · Jul 5, 2023 · Jul 5, 2023 · Jul 5, 2023
diff --git a/S1_NRB/metadata/extract.py b/S1_NRB/metadata/extract.py
@@ -265,6 +265,61 @@ def convert(obj, type):
         return out
 
 
+def calc_enl(tif, block_size=25, return_arr=False):
+    """
+    Calculate the equivalent number of looks (ENL) for a linear-scaled backscatter measurement GeoTIFF file of the
+    product scene.
+
+    Parameters
+    ----------
+    tif: str
+        The path to a linear-scaled backscatter measurement GeoTIFF file of the product scene.
+    block_size: int, optional
+        The block size to use for the calculation. Default is 25, which means that ENL will be calculated for 25x25
+        pixel blocks.
+    return_arr: bool, optional
+        If True, the calculated ENL array is returned. Default is False.
+
+    Returns
+    -------
+    float or numpy.ndarray
+        The median ENL value or array of ENL values if `return_enl_arr` is True.
+
+    References
+    ----------
+    .. [1]  S. N. Anfinsen, A. P. Doulgeris and T. Eltoft,
+            "Estimation of the Equivalent Number of Looks in Polarimetric Synthetic Aperture Radar Imagery,"
+            in IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 11, pp. 3795-3809, Nov. 2009,
+            doi: 10.1109/TGRS.2009.2019269.
+    """
+    with Raster(tif) as ras:
+        arr = ras.array()
+
+    arr = np.where(np.isinf(arr), np.nan, arr)
+
+    num_blocks_rows = arr.shape[0] // block_size
+    num_blocks_cols = arr.shape[1] // block_size
+    if num_blocks_rows == 0 or num_blocks_cols == 0:
+        raise ValueError("Block size is too large for the input data dimensions.")
+
+    blocks = arr[:num_blocks_rows * block_size,
+             :num_blocks_cols * block_size].reshape(
+        num_blocks_rows, block_size, num_blocks_cols, block_size
+    )
+
+    _mean = np.nanmean(blocks, axis=(1, 3))
+    _std = np.nanstd(blocks, axis=(1, 3))
+    enl = np.divide(_mean ** 2, _std ** 2, out=np.full_like(_mean, fill_value=np.nan), where=_std != 0)
+
+    out_arr = np.zeros((num_blocks_rows, num_blocks_cols))
+    out_arr[:num_blocks_rows, :num_blocks_cols] = enl
+
+    if return_arr:
+        return out_arr
+    else:
+        return np.nanmedian(out_arr)
+
+
 def calc_performance_estimates(files):
     """
     Calculates the performance estimates specified in CARD4L NRB 1.6.9 for all noise power images if available.
@@ -496,7 +551,7 @@ def meta_dict(config, target, src_ids, rtc_dir, proc_time, start, stop, compress
     sid0 = src_sid[list(src_sid.keys())[0]]  # first key/first file; used to extract some common metadata
     swath_id = re.search('_(IW|EW|S[1-6])_', os.path.basename(sid0.file)).group().replace('_', '')
 
-    ref_tif = finder(target, ['[hv]{2}-[gs]-lin.tif$'], regex=True)[0]
+    ref_tif = finder(target, ['(vh|hv)-[gs]-lin.tif$'], regex=True)[0]
     np_tifs = finder(target, ['-np-[hv]{2}.tif$'], regex=True)
     ei_tifs = finder(target, ['-ei.tif$'], regex=True)
     if len(ei_tifs) > 0:
@@ -575,6 +630,7 @@ def meta_dict(config, target, src_ids, rtc_dir, proc_time, start, stop, compress
     meta['prod']['demAccess'] = dem_access
     meta['prod']['doi'] = config['meta']['doi']
     meta['prod']['ellipsoidalHeight'] = None
+    meta['prod']['equivalentNumberLooks'] = np.round(calc_enl(tif=ref_tif), 2)
     meta['prod']['fileBitsPerSample'] = '32'
     meta['prod']['fileByteOrder'] = 'little-endian'
     meta['prod']['fileDataType'] = 'float'

diff --git a/S1_NRB/metadata/stac.py b/S1_NRB/metadata/stac.py
@@ -72,7 +72,8 @@ def product_json(meta, target, tifs, exist_ok=False):
                   polarizations=[Polarization[pol] for pol in meta['common']['polarisationChannels']],
                   product_type=meta['prod']['productName-short'],
                   looks_range=int(meta['prod']['rangeNumberOfLooks']),
-                  looks_azimuth=int(meta['prod']['azimuthNumberOfLooks']))
+                  looks_azimuth=int(meta['prod']['azimuthNumberOfLooks']),
+                  looks_equivalent_number=meta['prod']['equivalentNumberLooks'])
 
     sat_ext.apply(orbit_state=OrbitState[meta['common']['orbitDirection'].upper()],
                   relative_orbit=meta['common']['orbitNumbers_rel']['stop'],

diff --git a/S1_NRB/metadata/xml.py b/S1_NRB/metadata/xml.py
@@ -404,6 +404,8 @@ def product_xml(meta, target, tifs, nsmap, exist_ok=False):
     azimuthNumberOfLooks.text = str(meta['prod']['azimuthNumberOfLooks'])
     rangeNumberOfLooks = etree.SubElement(earthObservationMetaData, _nsc('s1-nrb:rangeNumberOfLooks', nsmap))
     rangeNumberOfLooks.text = str(meta['prod']['rangeNumberOfLooks'])
+    equivalentNumberLooks = etree.SubElement(earthObservationMetaData, _nsc('s1-nrb:equivalentNumberOfLooks', nsmap))
+    equivalentNumberLooks.text = str(meta['prod']['equivalentNumberLooks'])
     refDoc = etree.SubElement(earthObservationMetaData, _nsc('s1-nrb:refDoc', nsmap),
                               attrib={'name': meta['prod']['productName'],
                                       'version': meta['prod']['card4l-version'],

diff --git a/docs/_assets/enl_example_tile.png b/docs/_assets/enl_example_tile.png
diff --git a/docs/_assets/enl_grd_comparison_aois.png b/docs/_assets/enl_grd_comparison_aois.png
diff --git a/docs/_assets/enl_grd_comparison_scatter.png b/docs/_assets/enl_grd_comparison_scatter.png
diff --git a/docs/api.rst b/docs/api.rst
@@ -165,8 +165,10 @@ Extraction
     .. autosummary::
         :nosignatures:
 
+        calc_enl
         calc_geolocation_accuracy
         calc_performance_estimates
+        copy_src_meta
         etree_from_sid
         extract_pslr_islr
         find_in_annotation

diff --git a/docs/general/enl.rst b/docs/general/enl.rst
@@ -0,0 +1,75 @@
+Equivalent Number of Looks (ENL)
+================================
+
+The Equivalent Number of Looks (ENL) describes the degree of averaging applied to SAR measurements during data formation
+and postprocessing and is an indirect measure of speckle reduction (e.g., due to multilooking or speckle filtering).
+
+In case of linear scaled backscatter data, ENL can be calculated as:
+
+.. math::
+    ENL = \frac{\mu^2}{\sigma^2}
+
+where :math:`\mu` is the mean and :math:`\sigma` is the standard deviation of the image. ([1], section A1.1.7)
+
+The ENL value stored in the metadata of each S1-NRB product is calculated as suggested in [2], where ENL is first
+calculated for small pixel windows over the cross-polarized backscatter image and afterwards the median value of
+the distribution is selected.
+
+Calculate ENL per image
+-----------------------
+While only the median value is currently stored in the metadata of each S1-NRB product, it is possible to calculate ENL
+as described above for entire images using the function :func:`S1_NRB.metadata.extract.calc_enl`. The following code
+example shows how to calculate ENL for 25x25 pixel windows and return the result as a numpy array. The visualization of
+the resulting array is shown in Figure 1.
+
+.. code-block:: python
+
+    from S1_NRB.metadata.extract import calc_enl
+
+    tif = "s1a-iw-nrb-20220721t051225-044194-05465e-33tuf-vh-s-lin.tif"
+    enl_arr = calc_enl(tif=tif, block_size=25, return_arr=True)
+
+
+.. figure:: ../_assets/enl_example_tile.png
+    :width: 85 %
+    :align: center
+    :alt: Figure 1: Visualized ENL array for a S1-NRB product.
+
+    Figure 1: Visualized ENL array for a S1-NRB product processed from a Sentinel-1A SLC scene in IW mode for MGRS tile 33TUF
+    (coastline between Rome and Naples, Italy).
+
+Comparison between GRDH and NRB
+-------------------------------
+[1] provides estimates of ENL for different Sentinel-1 products (average over all swaths), e.g. ENL of 4.4 for GRDH in
+IW mode, and a description of the estimation process in section D1. The following shows a simple comparison between the
+GRDH product:
+
+`S1A_IW_GRDH_1SDV_20220721T051222_20220721T051247_044194_05465E_5807`
+
+and a S1-NRB product derived from the equivalent SLC product and processed for MGRS tile 33TUF:
+
+`S1A_IW_SLC__1SDV_20220721T051221_20220721T051249_044194_05465E_BACD`
+
+ENL was calculated for a selection of homogeneous forest areas, which are highlighted in Figure 2. The green outline
+traces the north-western corner of MGRS tile 33TUF (see Fig. 1). The resulting scatter plot (Figure 3) shows
+consistently higher ENL values for the GRDH product (Avg. ENL: 4.81) in comparison to the S1-NRB product (Avg. ENL: 4.59).
+
+.. figure:: ../_assets/enl_grd_comparison_aois.png
+    :width: 75 %
+    :align: center
+    :alt: Figure 2: Selection of homogeneous forest areas for ENL comparison between GRDH and NRB.
+
+    Figure 2: Selection of homogeneous forest areas for ENL comparison between GRDH and NRB. Green outline: North-western
+    corner of MGRS tile 33TUF; Background image: VH backscatter of the GRDH product.
+
+.. figure:: ../_assets/enl_grd_comparison_scatter.png
+    :width: 75 %
+    :align: center
+    :alt: Figure 3: Scatter plot comparing ENL values between GRDH and NRB, calculated for selected areas (see Fig. 2).
+
+    Figure 3: Scatter plot comparing ENL values between GRDH and NRB, calculated for selected areas (see Fig. 2).
+
+References
+----------
+* [1] CLS (2016): Sentinel-1 Product Definition. Version 2.7. URL: https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/document-library/-/asset_publisher/1dO7RF5fJMbd/content/sentinel-1-product-definition.
+* [2] S. N. Anfinsen, A. P. Doulgeris and T. Eltoft, "Estimation of the Equivalent Number of Looks in Polarimetric Synthetic Aperture Radar Imagery," in IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 11, pp. 3795-3809, Nov. 2009, doi: `10.1109/TGRS.2009.2019269 <https://doi.org/10.1109/TGRS.2009.2019269>`_.
diff --git a/docs/general/index.rst b/docs/general/index.rst
@@ -7,5 +7,6 @@ General Topics
     installation.rst
     usage.rst
     S1_NRB.rst
-    geoaccuracy.rst
     folderstructure.rst
+    geoaccuracy.rst
+    enl.rst