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,69 @@
+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`:
+
+.. code-block:: python
+
+    from S1_NRB.metadata.extract import calc_enl
+
+    enl_arr = calc_enl(tif="s1a-iw-nrb-20220721t051225-044194-05465e-33tuf-vh-s-lin.tif", block_size=25, return_arr=True)
+
+
+Example calculated from sigma0 VH backscatter of a S1-NRB product, which has been derived from a Sentinel-1 SLC (IW) scene:
+
+.. figure:: ../_assets/enl_example_tile.png
+    :width: 85 %
+    :align: center
+    :alt: ENL calculation for an example scene
+
+
+Comparison between GRDH and NRB
+-------------------------------
+[1] provides estimates of ENL for different Sentinel-1 products, e.g. 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:
+
+`S1A_IW_SLC__1SDV_20220721T051221_20220721T051249_044194_05465E_BACD`
+
+ENL was calculated for a selection of homogeneous forest areas, as can be seen in the following figure. The
+green outline traces the north-western corner of MGRS tile 33TUF (see also figure above).
+
+.. figure:: ../_assets/enl_grd_comparison_aois.png
+    :width: 75 %
+    :align: center
+    :alt: Selected areas for ENL comparison between GRDH and NRB
+
+
+The resulting scatter plot shows slightly higher ENL values for the GRDH product:
+
+.. figure:: ../_assets/enl_grd_comparison_scatter.png
+    :width: 75 %
+    :align: center
+    :alt: ENL comparison between GRDH and NRB
+
+
+References
+----------
+* [1] `Sentinel-1 Product Definition <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