From be20d258d678e4f1e8435856a9915bc4990dcfd8 Mon Sep 17 00:00:00 2001 From: telespazio-tim Date: Tue, 23 Apr 2024 14:15:44 +0200 Subject: [PATCH] Fix calculation of sen2cor input image ROI center --- sen2like/release-notes.md | 6 ++ sen2like/sen2like/core/product_preparation.py | 2 +- .../core/sen2cor_client/sen2cor_client.py | 89 ++++++++++++++++--- sen2like/sen2like/grids/grids.py | 10 --- sen2like/sen2like/grids/mgrs_framing.py | 30 ------- 5 files changed, 82 insertions(+), 55 deletions(-) diff --git a/sen2like/release-notes.md b/sen2like/release-notes.md index 38cfaad..73ab831 100644 --- a/sen2like/release-notes.md +++ b/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 diff --git a/sen2like/sen2like/core/product_preparation.py b/sen2like/sen2like/core/product_preparation.py index b169b21..467f1db 100644 --- a/sen2like/sen2like/core/product_preparation.py +++ b/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) diff --git a/sen2like/sen2like/core/sen2cor_client/sen2cor_client.py b/sen2like/sen2like/core/sen2cor_client/sen2cor_client.py index df6509b..0166719 100644 --- a/sen2like/sen2like/core/sen2cor_client/sen2cor_client.py +++ b/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 diff --git a/sen2like/sen2like/grids/grids.py b/sen2like/sen2like/grids/grids.py index 2b8b217..731c52b 100644 --- a/sen2like/sen2like/grids/grids.py +++ b/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. diff --git a/sen2like/sen2like/grids/mgrs_framing.py b/sen2like/sen2like/grids/mgrs_framing.py index 8e57f5c..306cdd1 100644 --- a/sen2like/sen2like/grids/mgrs_framing.py +++ b/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