Merge branch 'develop' into 'master' to release 4.4.3

Author: telespazio-tim

sen2like/release-notes.md

@@ -1,5 +1,11 @@
 # Sen2Like Release Notes
 
+## v4.4.3
+
+### Fix
+
+* Fix calculation of Sen2cor region of interest when target UTM differs from the input Landsat product UTM.
+
 ## v4.4.2
 
 ### Release of Sen2Cor3 documentation and software version 3.01.00

sen2like/sen2like/core/product_preparation.py

@@ -318,7 +318,7 @@ def _extract_product_files(self, product: S2L_Product):
         if (
             scl_dir
             and (not product.context.use_sen2cor)
-            and product_reader.data_type != "Level-2A"
+            and product_reader.data_type != "Level-2A" # Sentinel
         ):
             product_reader.scene_classif_band = self._get_scl_map(scl_dir, product)
 

sen2like/sen2like/core/sen2cor_client/sen2cor_client.py

@@ -16,17 +16,41 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-
+import affine
 import logging
 import os
 import subprocess
 import lxml.etree as ET
+
+from osgeo import osr
+import mgrs
+
 from core import S2L_config
-from grids.mgrs_framing import pixel_center
+from core.image_file import S2L_ImageFile
+
 
 logger = logging.getLogger("Sen2Like")
 
 
+def get_mgrs_center(tilecode: str, utm=False) -> tuple:
+    """Get MGRS tile center coordinates in native tile UTM or WGS84
+
+    Args:
+        tilecode (str): tile to get center coords
+        utm (bool, optional): Flag for output SRS . Defaults to False.
+
+    Returns:
+        tuple: (lat,long) coords if utm=False, else (utm, N/S, easting, northing)
+    """
+    if tilecode.startswith('T'):
+        tilecode = tilecode[1:]
+    centercode = tilecode + '5490045100'
+    m = mgrs.MGRS()
+    if utm:
+        return m.MGRSToUTM(centercode)
+    return m.toLatLon(centercode)
+
+
 class Sen2corClient:
 
     gipp_template_file = os.path.join(
@@ -117,26 +141,63 @@ def _write_gipp(self, product):
         # ref_band = None is considered as S2 product format (S2A, S2B, S2P prisma)
         ref_band = self.roi_ref_band.get(product.mtl.mission, None)
 
-        if ref_band is None:
-            logger.debug("For sentinel, sen2cor don't use ROI")
-            ET.ElementTree(root).write(gipp_path, encoding='utf-8', xml_declaration=True)
-            return gipp_path
+        if ref_band:
+            # Compute ROI center for landsat and fill template with result
+            ref_band_file = product.get_band_file(ref_band)
 
-        ref_band_file = product.get_band_file(ref_band)
+            y, x = self._pixel_center(ref_band_file)
 
-        y, x = pixel_center(ref_band_file, self.out_mgrs)
-        logger.debug('Pixel center : (%s, %s)', y, x)
+            logger.debug('Pixel center : (%s, %s)', y, x)
 
-        row0 = root.find('Common_Section/Region_Of_Interest/row0')
-        row0.text = str(y)
-        col0 = root.find('Common_Section/Region_Of_Interest/col0')
-        col0.text = str(x)
+            row0 = root.find('Common_Section/Region_Of_Interest/row0')
+            row0.text = str(y)
+            col0 = root.find('Common_Section/Region_Of_Interest/col0')
+            col0.text = str(x)
 
-        ET.ElementTree(root).write(gipp_path, encoding='utf-8', xml_declaration=True)
+            ET.ElementTree(root).write(gipp_path, encoding='utf-8', xml_declaration=True)
+
+        else:
+            logger.debug("For sentinel, sen2cor don't use ROI")
 
         logger.info('GIPP L2A : %s', gipp_path)
+
+        ET.ElementTree(root).write(gipp_path, encoding='utf-8', xml_declaration=True)
+
         return gipp_path
 
+    def _pixel_center(self, image: S2L_ImageFile):
+        """Get mgrs tile center coordinates from tile code in image coordinates
+
+        Args:
+            image (S2L_ImageFile): image to get srs from
+
+        Returns:
+            tuple: (y,x) image coordinates
+        """
+
+        lat, lon = get_mgrs_center(self.out_mgrs) # pylint: disable=W0632
+
+        # Transform src SRS
+        wgs84_srs = osr.SpatialReference()
+        wgs84_srs.ImportFromEPSG(4326)
+
+        # Transform dst SRS
+        image_srs = osr.SpatialReference(wkt=image.projection)
+
+        # convert MGRS center coordinates from lat lon to image EPSG coordinates (UTM)
+        transformation = osr.CoordinateTransformation(wgs84_srs, image_srs)
+        easting, northing, _ = transformation.TransformPoint(lat, lon)
+
+        # northing = y = latitude, easting = x = longitude
+        tr = affine.Affine(
+            image.yRes, 0, image.yMax,
+            0, image.xRes, image.xMin
+        )
+
+        # compute y,x in image coordinates
+        y, x = (northing, easting) * (~ tr)
+        return int(y), int(x)
+
 
 class Sen2corError(Exception):
     pass

sen2like/sen2like/grids/grids.py

@@ -22,7 +22,6 @@
 import sqlite3
 from os.path import abspath, dirname
 
-import mgrs
 import pandas as pd
 from osgeo import ogr, osr
 from shapely.wkt import loads
@@ -59,15 +58,6 @@ def getROIfromMGRS(self, tilecode):
         # return as dict
         return roi.to_dict(orient='list')
 
-    def get_mgrs_center(self, tilecode, utm=False):
-        if tilecode.startswith('T'):
-            tilecode = tilecode[1:]
-        centercode = tilecode + '5490045100'
-        m = mgrs.MGRS()
-        if utm:
-            return m.MGRSToUTM(centercode)
-        return m.toLatLon(centercode)
-
     # Don't know why but with this method the SRS code is not added in the geojson
     # if the GDAL_DATA variable is set in the environment. However we need
     # this GDAL_DATA variable for other methods. So we do not use this method.

sen2like/sen2like/grids/mgrs_framing.py

@@ -24,7 +24,6 @@
 from math import ceil
 from typing import NamedTuple
 
-import affine
 import numpy as np
 import pandas as pd
 from numpy.typing import NDArray
@@ -129,35 +128,6 @@ def resample(image: S2L_ImageFile, res: int, filepath_out: str) -> S2L_ImageFile
     return image.duplicate(filepath_out, array=data, res=res)
 
 
-def pixel_center(image: S2L_ImageFile, tile_code: str):
-    """Get mgrs tile center coordinates from tile code in image SRS
-
-    Args:
-        image (S2L_ImageFile): image to get srs from
-        tile_code (str): MGRS tile code
-
-    Returns:
-        tuple: y/x coordinates
-    """
-    converter = grids.GridsConverter()
-    utm, orientation, easting, northing = converter.get_mgrs_center(tile_code, utm=True)
-    # UTM South vs. UTM North ?
-    tile_srs = osr.SpatialReference()
-    tile_srs.ImportFromEPSG(int('32' + ('6' if orientation == 'N' else '7') + str(utm)))
-    image_srs = osr.SpatialReference(wkt=image.projection)
-    if not tile_srs.IsSame(image_srs):
-        transformation = osr.CoordinateTransformation(tile_srs, image_srs)
-        northing, easting = transformation.TransformPoint((northing, easting))
-
-    # northing = y = latitude, easting = x = longitude
-    tr = affine.Affine(
-        image.yRes, 0, image.yMax,
-        0, image.xRes, image.xMin
-    )
-    y, x = (northing, easting) * (~ tr)
-    return int(y), int(x)
-
-
 def _reproject(filepath: str, dir_out: str, ds_src: gdal.Dataset, x_res: int, y_res: int, target_srs: osr.SpatialReference, order: int) -> ReprojectResult:
     """Reproject given dataset in the target srs in the given resolution