Merge pull request #40 from SAR-ARD/geolocation_accuracy

Geolocation accuracy

S1_NRB/metadata/extract.py

@@ -1,6 +1,6 @@
 import os
 import re
-from math import isclose
+import math
 from lxml import etree
 from datetime import datetime
 import numpy as np
@@ -12,7 +12,7 @@
 from spatialist.raster import rasterize
 from osgeo import gdal
 import S1_NRB
-from S1_NRB.metadata.mapping import NRB_PATTERN, ITEM_MAP, RES_MAP, ORB_MAP, DEM_MAP
+from S1_NRB.metadata.mapping import NRB_PATTERN, ITEM_MAP, RES_MAP, ORB_MAP, DEM_MAP, SLC_ACC_MAP
 
 gdal.UseExceptions()
 
@@ -100,10 +100,10 @@ def vec_from_srccoords(coord_list):
     """
     if len(coord_list) == 2:
         # determine joined border between footprints
-        if isclose(coord_list[0][0][0], coord_list[1][3][0], abs_tol=0.1):
+        if math.isclose(coord_list[0][0][0], coord_list[1][3][0], abs_tol=0.1):
             c1 = coord_list[1]
             c2 = coord_list[0]
-        elif isclose(coord_list[1][0][0], coord_list[0][3][0], abs_tol=0.1):
+        elif math.isclose(coord_list[1][0][0], coord_list[0][3][0], abs_tol=0.1):
             c1 = coord_list[0]
             c2 = coord_list[1]
         else:
@@ -427,6 +427,60 @@ def _get_block_offset(band):
     return headers_size
 
 
+def calc_geolocation_accuracy(swath_identifier, ei_tif, dem_type, etad):
+    """
+    Calculates the radial root mean square error, which is a target requirement of the CARD4L NRB specification
+    (Item 4.3). For more information see: https://s1-nrb.readthedocs.io/en/latest/general/geoaccuracy.html
+    
+    Parameters
+    ----------
+    swath_identifier: str
+        Swath identifier dependent on acquisition mode.
+    ei_tif: str
+        Path to the annotation GeoTIFF layer 'Ellipsoidal Incident Angle' of the current product.
+    dem_type: str
+        The DEM type used for processing.
+    etad: bool
+        Was the ETAD correction applied?
+    
+    Returns
+    -------
+    rmse_planar: float
+        The calculated rRMSE value rounded to two decimal places.
+    """
+    if 'copernicus' not in dem_type.lower():
+        return None
+
+    if etad:
+        # https://sentinel.esa.int/nl/web/sentinel/missions/sentinel-1/data-products/etad-dataset
+        slc_acc = {'ALE': {'rg': 0,
+                           'az': 0},
+                   '1sigma': {'rg': 0.2,
+                              'az': 0.1}}
+    else:
+        swath_id = 'SM' if re.search('S[1-6]', swath_identifier) is not None else swath_identifier
+        slc_acc = SLC_ACC_MAP[swath_id]
+
+    # min/max ellipsoidal incidence angle
+    with Raster(ei_tif) as ras:
+        stats = ras.allstats(approximate=False)
+        ei_min = stats[0]['min']
+
+    # COP-DEM global mean accuracy (LE68) based on LE90 under assumption of gaussian distribution:
+    copdem_glob_1sigma_le68 = 1.56
+    rmse_dem_planar = copdem_glob_1sigma_le68 / math.tan(math.radians(ei_min))
+
+    # Calculation of RMSE_planar
+    rmse_rg = math.sqrt(slc_acc['ALE']['rg'] ** 2 + slc_acc['1sigma']['rg'] ** 2)
+    rmse_az = math.sqrt(slc_acc['ALE']['az'] ** 2 + slc_acc['1sigma']['az'] ** 2)
+
+    rmse_planar = math.sqrt((rmse_rg / math.sin(math.radians(ei_min))) ** 2 +
+                            rmse_az ** 2 +
+                            rmse_dem_planar ** 2)
+
+    return round(rmse_planar, 2)
+
+
 def meta_dict(config, target, src_ids, snap_datasets, proc_time, start, stop, compression):
     """
     Creates a dictionary containing metadata for a product scene, as well as its source scenes. The dictionary can then
@@ -468,9 +522,11 @@ def meta_dict(config, target, src_ids, snap_datasets, proc_time, start, stop, co
         src_sid[uid] = sid
         src_xml[uid] = etree_from_sid(sid=sid)
     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}-g-lin.tif$'], regex=True)[0]
     np_tifs = finder(target, ['-np-[hv]{2}.tif$'], regex=True)
+    ei_tif = finder(target, ['-ei.tif$'], regex=True)[0]
 
     product_id = os.path.basename(target)
     prod_meta = get_prod_meta(product_id=product_id, tif=ref_tif, src_ids=src_ids,
@@ -490,6 +546,9 @@ def meta_dict(config, target, src_ids, snap_datasets, proc_time, start, stop, co
     tups = [(ITEM_MAP[key]['suffix'], ITEM_MAP[key]['z_error']) for key in ITEM_MAP.keys()]
     z_err_dict = dict(tups)
 
+    geocorr_acc = calc_geolocation_accuracy(swath_identifier=swath_id, ei_tif=ei_tif, dem_type=dem_type,
+                                            etad=config['etad'])
+
     # Common metadata (sorted alphabetically)
     meta['common']['antennaLookDirection'] = 'RIGHT'
     meta['common']['constellation'] = 'sentinel-1'
@@ -514,8 +573,7 @@ def meta_dict(config, target, src_ids, snap_datasets, proc_time, start, stop, co
     meta['common']['radarBand'] = 'C'
     meta['common']['radarCenterFreq'] = 5405000000
     meta['common']['sensorType'] = 'RADAR'
-    meta['common']['swathIdentifier'] = re.search('_(IW|EW|S[1-6])_',
-                                                  os.path.basename(sid0.file)).group().replace('_', '')
+    meta['common']['swathIdentifier'] = swath_id
     meta['common']['wrsLongitudeGrid'] = str(sid0.meta['orbitNumbers_rel']['start'])
 
     # Product metadata (sorted alphabetically)
@@ -551,8 +609,9 @@ def meta_dict(config, target, src_ids, snap_datasets, proc_time, start, stop, co
     meta['prod']['geoCorrAccuracyEasternSTDev'] = None
     meta['prod']['geoCorrAccuracyNorthernBias'] = None
     meta['prod']['geoCorrAccuracyNorthernSTDev'] = None
-    meta['prod']['geoCorrAccuracy_rRMSE'] = None
-    meta['prod']['geoCorrAccuracyReference'] = None
+    meta['prod']['geoCorrAccuracy_rRMSE'] = geocorr_acc
+    meta['prod']['geoCorrAccuracyReference'] = 'https://s1-nrb.readthedocs.io/en/{}/general/geoaccuracy.html' \
+                                               ''.format(S1_NRB.__version__)
     meta['prod']['geoCorrAccuracyType'] = 'slant-range'
     meta['prod']['geoCorrAlgorithm'] = 'https://sentinel.esa.int/documents/247904/1653442/' \
                                        'Guide-to-Sentinel-1-Geocoding.pdf'

S1_NRB/metadata/mapping.py

@@ -150,3 +150,17 @@
                                   'unit': 'dB',
                                   'role': 'noise-power',
                                   'title': 'Noise Power'}}
+
+# https://sentinel.esa.int/documents/247904/1653442/Guide-to-Sentinel-1-Geocoding.pdf
+SLC_ACC_MAP = {'SM': {'ALE': {'rg': -3.02,
+                              'az': 2.02},
+                      '1sigma': {'rg': 0.26,
+                                 'az': 0.41}},
+               'IW': {'ALE': {'rg': -2.99,
+                              'az': 2.03},
+                      '1sigma': {'rg': 0.22,
+                                 'az': 0.53}},
+               'EW': {'ALE': {'rg': -3.41,
+                              'az': 1.37},
+                      '1sigma': {'rg': 0.7,
+                                 'az': 2.27}}}

S1_NRB/metadata/xml.py

@@ -433,10 +433,10 @@ def product_xml(meta, target, tifs, nsmap, exist_ok=False):
     geoCorrAccuracyEasternBias.text = meta['prod']['geoCorrAccuracyEasternBias']
     geoCorrAccuracy_rRMSE = etree.SubElement(earthObservationMetaData,
                                              _nsc('s1-nrb:geoCorrAccuracy_rRMSE', nsmap), attrib={'uom': 'm'})
-    geoCorrAccuracy_rRMSE.text = meta['prod']['geoCorrAccuracy_rRMSE']
+    geoCorrAccuracy_rRMSE.text = str(meta['prod']['geoCorrAccuracy_rRMSE'])
+    geoacc_ref = meta['prod']['geoCorrAccuracyReference']
     geoCorrAccuracyReference = etree.SubElement(earthObservationMetaData, _nsc('s1-nrb:geoCorrAccuracyReference', nsmap),
-                                                attrib={_nsc('xlink:href', nsmap): meta['prod'][
-                                                    'geoCorrAccuracyReference']})
+                                                attrib={_nsc('xlink:href', nsmap): geoacc_ref})
     numLines = etree.SubElement(earthObservationMetaData, _nsc('s1-nrb:numLines', nsmap))
     numLines.text = meta['prod']['numLines']
     numPixelsPerLine = etree.SubElement(earthObservationMetaData, _nsc('s1-nrb:numPixelsPerLine', nsmap))

docs/api.rst

@@ -123,6 +123,7 @@ Extraction
     .. autosummary::
         :nosignatures:
 
+        calc_geolocation_accuracy
         calc_performance_estimates
         convert_coordinates
         etree_from_sid

docs/general/geoaccuracy.rst

@@ -0,0 +1,48 @@
+Geolocation Accuracy
+====================
+
+Item 4.3 of the CARD4L NRB specification requires, as minimum, an estimate of the absolute location error (ALE) “as
+bias and standard deviation, provided in slant range/azimuth, or Northing/Easting” [1]. As desired target the accuracy
+is less or equal 0.1 pixels radial root mean square error (rRMSE), which can be defined as:
+
+.. math::
+   RMSE_{planar} = \sqrt{RMSE_{SLC,Az}^2 + (\frac{RMSE_{SLC,Rg}}{sin(\theta_{i,min})})^2 + RMSE_{DEM,planar}^2 + RMSE_{proc}^2}
+
+The error induced by the DEM can be described as:
+
+.. math::
+   RMSE_{DEM,planar} = \frac{\sigma_{DEM}}{tan(\theta_{i,min})}
+
+where
+
+:math:`\theta_{i,min}` = The minimum possible angle of incidence
+
+:math:`RMSE_{SLC,Az/Rg}` = Error induced by SLC source data in azimuth/range
+
+:math:`RMSE_{DEM,planar}` = Error induced by DEM inaccuracy
+
+:math:`RMSE_{proc}` = Error induced by other processing steps
+
+:math:`\sigma_{DEM}` = DEM accuracy at :math:`1\sigma` (LE68)
+
+
+Limitations
+-----------
+Currently, the following simplifications need to be considered for the calculation of rRMSE values found in the metadata
+of each S1-NRB product:
+
+* Processing induced errors (:math:`RMSE_{proc}`) and the error term related to DEM interpolation are not further considered and assumed to be 0.
+* The DEM accuracy (:math:`\sigma_{DEM}`) is estimated on the global mean accuracy LE90 reported for the COP-DEM [2] under the assumption of gaussian distribution:
+
+    * Global: LE90 = 2.57; LE68 :math:`\approx` 1.56
+
+* rRMSE is only calculated if a COP-DEM was used for processing, otherwise the value is set to ``None``
+
+Development Status
+------------------
+The development status is tracked and discussed in the following Github issue: https://github.com/SAR-ARD/S1_NRB/issues/33
+
+References
+-----------
+* [1] https://ceos.org/ard/files/PFS/NRB/v5.5/CARD4L-PFS_NRB_v5.5.pdf
+* [2] https://spacedata.copernicus.eu/documents/20126/0/GEO1988-CopernicusDEM-SPE-002_ProductHandbook_I1.00.pdf/082dd479-f908-bf42-51bf-4c0053129f7c?t=1586526993604

docs/general/index.rst

@@ -7,3 +7,4 @@ General Topics
    installation.rst
    usage.rst
    S1_NRB.rst
+   geoaccuracy.rst